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ABSTRACT 

The  effect  of  added  resistance  due  to  steering  on  a 
high-speed  containership  is  propulsion  reduction. 
Limitation  on  the  propulsion  losses  can  be  achieved  by  a 
properly  designed  controller,  which  minimizes  the  rudder 
activity  as  well  as  providing  desired  overseas  heading. 

A  computer  program,  simulating  a  cascade  configuration 
of  the  SL-7  high-speed  containership  along  with  a  specific 
controller  was  coupled  to  a  function  minimization  subroutine 
as  well  as  a  sea  state  generator  subroutine  in  order  to 
minimize  a  performance  criterion. 

The  entire  model  was  tested  for  a  fixed  speed,  several 
encounter  frequencies,  several  encounter  angles  in  calm 
waters  and  in  a  seaway  as  well. 
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I.  INTRODUCTION 

The  1973  oil  crisis, where  the  price  of  a  barrel  of  oil 
jumped  higher  than  210  percent  and  affected  every  aspect  of 
the  world  economy , directly  affected  the  shipping  transporta- 
tion industries .Whereas  the  price  of  oil  was  formerly  of 
modest  importance,  it  became  a  prime  concern. Extensive  asso- 
ciation to  the  emerged  problem  has  to  do  with  the  design  of 
autopilots  for  ships,  which  can  minimize  propulsive  losses 
caused  by  added  resistance  due  to  steering. 

One  motivation  for  the  design  of  an  optimum  steering 
controller  was  the  work  done  by  Nomoto  and  Motoyama  who 
claimed  that  reduction  in  propulsion  loss  could  amount  to  a 
1  percent  savings  in  fuel  consumption. 

For  most  commercial  ships  a  1  to  2  percent  saving  in 
fuel  costs , justifies  the  expense  of  fitting  an  autopilot 
which  has  the  capability  of  producing  this  savings  [Ref.  1] . 

Chapter  2  addresses  what  type  of  computer  model  can  be 
used  to  represent  the  ship. Several  of  these  models  are 
investigated, such  as  simulation  from  the  equations  of  motion 
and  the  Nomoto  third  order  which  was  developed  from  the 
equation  of  motion. 

Chapter  3  addresses  the  problem  of  selecting  an  adequate 
cost  function  to  represent  the  added  resistance  due  to 
steering. 

The  problem  of  finding  the  best  controller  design  to 
provide  a  minimum  value  of  added  resistance  due  to  steering, 
for  regular  and  irregular  seas, is  studied  in  Chapters  4  and 
5  respectively  by  using  as  a  model  of  the  ship  the  equations 
of  motion  and  a  function  minimization  subroutine. 

Chapter  6  indicates  how  the  controller  parameters  can  be 
adjusted  in  any  encounter  environmental  condition  by  means 
of  an  adaptive  control. 

11 


Chapter  7  presents  the  conclusions  from  the  experimental 
work. 
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II.  DESCRIPTION  OF  COMPUTER  MODELS 

The  most  accurate  model  which  can  represent  the  ship/ 
steering  dynamics  is  the  model  which  is  based  upon  the  equa- 
tions of  motion  as  defined  by  series  expansion  including  all 
terms  (both  linear  and  nonlinear) . 

Using  experimentally  measured  hydrodynamic  coefficients, 
for  the  SL-7  high  speed  containership ,  a  computer  program 
was  developed  in  order  to  provide  a  computer  simulation  for 
the  ship. Figure  2.1  shows  the  block  diagram  and  the  computer 
program  can  be  found  in  Appendix  A  [Ref.  2,3]. 
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Figure  2.1    Block  Diagram  of  Ship  and  Control  System 
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Figure  2.2    Determination  of  Nomoto  Third  Order  Model 

As  a  next  choice,  from  the  equations  of  motion  we  derive 
the  Nomoto  third  order  transfer  function.  Figure  2.2  shows 
the  block  diagram. 

Using  the  scheme  of  Figure  2.2  which  includes  the  func- 
tion minimization  subroutine  with  both  yaw  and  sway  equa- 
tions, with  the  linear  terms  only  we  obtain  the  appropriate 
coefficients  of  the  Nomoto  third  model.  (Including  in  the 
equations  of  motion  nonlinear  terms  we  can  see  that  the 
perturbations  were  small  enough). 

The  function  minimization  subroutine  used  was  BOXPLX, 
which  was  programed  by  R.  Hilleary.The  task  of  the  already 
mentioned  subroutine  is  to  find  the  minimum  of  any  function 
and  is  subjected  to  explicit  constraints  of  the  variables  or 
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implicit  constraints  on  functions   of  the  variables. In  addi- 
tion it  can  handle  a  maximum  of  25  variables  [Ref.  4]. 

The  extracted  results  were  very  close  to  the  results  we 
got  by  the  analytic  solution  and  are  tabulated  below  only 
for  16,  23  and  32  knots. 
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III.  STEERING  PERFORMANCE  CRITERION 

Studying  the  literature,  one  can  see  many  approaches  to 
the  problem  of  optimizing  an  automatic  ship  steering 
controller  for  maximum  reduction  in  fuel  consumption.  Since 
added  resistance  is  directly  related  to  both  rudder  activity 
and  yawing  motion  we  can  express  a  measure  of  this  added 
resistance  given  as  a  performance  criterion  with  the 
formula : 


(3.1) 


Where : 

•  Delta  (o   )=rudder  angle 

•  Psi  (fi)      =yaw  angle 

•  Lambda ( y^ )= weighting  factor 

This  formula  defines  an  approximate  drag  due  to  steering 
for  small  amplitude  oscillations  about  a  steady-state  pivot 
point  of  the  ship  during  yawing  at  the  natural  frequency  of 
the  closed-loop  ship/ steering  system. (About  0.05  rad/sec  for 
the  SL-7  containership) . It  is  also  convenient  for  shipboard 
use  because  yaw  error  and  rudder  angle  can  be  easily  meas- 
ured. 

During  this  study  the  values  for  the  weighting  factors 
are  taken  from  R.  E.  Reid's  work  for  the  high-speed  contain- 
ership SL-7  [Ref .  5] . 
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TABLE 

2 

Weighting 

Factors 

ts 

Lambda 

16, 

,796 

8, 

.128 

4, 

.2 

Ship's  Speed, Knots 

16 
23 
32 


In  Table  2  weighting  factors   for  the  operating  range  of 
the  ship  are  tabulated. 
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IV.  CONTROLLER  DESIGN  FOR  REGULAR  SEAS 

Our  goal  now,  is  to  estimate  the  system's  performance  in 
a  seaway.  To  accomplish  the  above  task  first  we  must  deter- 
mine a  suitable  representation  of  the  external  disturbances 
on  the  ship  by  the  sea.  This  can  be  done,  since  a  suffi- 
ciently accurate  computer  ship  model  and  a  steering  perform- 
ance criterion  have  been  defined. 

It  can  be  postulated  that  a  sufficiently  accurate  model 
of  the  seaway  itself,  is  a  representative  modeling  of  forces 
and  moments  exerted  on  the  ship  [Ref.  6,7,8].  An  explana- 
tion of  what  a  sea  comprises,  and  how  predicted  or  observed 
sea  states  can  be  analyzed  to  determine  the  forces  and 
motions  of  a  body  in  a  sea  was  studied  by  Michael  [Ref.  9]. 

In  this  chapter  we  will  use  as  a  seaway  representation 
the  regular  sea  model  in  which  the  forces  on  the  ship 
[Ref.  10],  have  the  form: 


W^'  R^'COs(ci)^i   +l?A 


(4.1) 


Where : 

•  Rc  =exciting  force 
•Cue  =encounter  frequency 

•  VVh  =wave  height 

•  ^i   = phase  angle 

Values  for  the  exciting  force  (  )  for  different 
encounter  angles  and  different  encounter  frequencies  as 
well,  were  taken  from  [Ref.  3].  Table  3  shows  the  corre- 
spondence between  wave  height  and  sea  state  which  were  used 
in  this  work. 
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The  controller   used  in   the  present   study  is   shown  in 
Appendix  B  and  its  form  is   repeated  here  in  Figure  4.1  . 

TABLE  3 
Sea  State  Vs  Wave  Height 


Sea  State 

(Beaufort  scale) 
1 
2 
4 
6 
7 


Wave  Height 

(in  feet) 

4.0 

4.5 

5.5 
10.0 
17.5 


— 

K2(1+T;^S)     (l+Ti^S) 
(I+T2S)     (I+T3G) 

Figure  4 . 1 


Controller  C 


Values  for   the  optimal  gains   for  the   above  controller 
and  values  for  the  cost  J  as  well,  are  shown  below  in  Tables 

4,5,6,7,8,9,10,11 

These  values  were  obtained  using: 

•  Sea  states  1-2-4-6-7  (Beaufort  scale). 

•  Encounter  angles  0° -  30° -  60° -  90° . 

•  Encounter  frequencies  0.2-0.4-0.6-0.75-1.5  rad  per  sec. 

•  Constant  speed  23  knots, 

A  careful  analysis   of  the  extracted  results   lead  us  to 
conclude : 
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•  The  maximum  deviation  of  controller  parameter  values 
occured  at  0.4  rad  per  sec  encounter  frequency  for  all 
used  encounter  angles  and  sea  states. 

•  For  the  same  encounter  f requency , f or  all  encounter 
angles  studied,  the  controller  parameter  values 
increasing  smoothly  as  the  sea  state  increases. 

•  For  all  encounter  frequencies  and  for  specific 
encounter  angles,  we  observe  that  the  cost  increases  at 
higher  sea  states. 

•  For  0.4  rad  per  sec  encounter  frequency  the  cost 
changes  rapidly  for  sea  state  6  and  7  as  we  move 
through  the  encounter  angles . 

•  For  the  same  encounter  angles  and  sea  states  the  cost 
decreases  at  higher  encounter  frequencies. 

•  For  all  encounter  angles  and  encounter  frequencies  the 
maximum  cost  occured  for  sea  state  7. 

•  The  cost  had  it's  maximum  value  for  encounter  angle  90° 
and  encounter  frequency  0.4  rad  per  sec. 

Furthermore  in  order  to  obtain  the  behavior  of  the 
rudder  and  yaw  motion  of  the  ship,  transient  response  plots 
were  obtained  for  controller  'C'  at  ship's  speed  23  knots, 
different  sea  states  and  encounter  angles  as  shown  in 
Figures  4.2  through  4.15 

These  plots  were  obtained  using  the  program  of  Appendix 
C. 

From  these  plots  it  is  verified  that  for  the  same 
encounter  frequencies  and  encounter  angles ,  rudder  and  yaw 
motion  increases  in  amplitude  as  the  sea  states  rises. 

Finally  Figures  4.16,4.17,4.18,4.19,4.20,  show  the 
results  of  an  experiment  in  which  we  were  changing  one 
controller  parameter  keeping  the  rest  fixed,  and  we  observed 
that  the  cost  does  not  change  significantly  in  the  vicinity 
of  the  actual  value. 
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Using  that  as  a  fact  we  can  postulate  that  for  a 
specific  sailing  mode  for  the  controller  parameter  values 
high  accuracy  is  not  required. 

Comparing  controller  C  with  controller  A  it  is  obvious 
that  the  parameter  values  of  controller  A  vary  over  wider 
range . 
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TABLE  4 

Controller  C  For  0°  Encounter  Wave  Angle 

Encounter  Frequency  0.2  Rads  Per  Sec 

Sea  State:            1  2 

Kl=             1.3999996  1.3999996 

Tl=            58.2299957  58.2299957 

T2=            24.0000000  24.0000000 

T3=            10.5000000  10.5000000 

T4=             3.0000000  3.0000000 

Cost  J=       0.6382236E-33  0 . 8077518E-33 

Encounter  Frequency  0.4  Rads  Per  Sec 

Sea  State:            1  2 

Kl=             1.3999996  1.3999996 

Tl=            58.2299957  58.2299957 

T2=            24.0000000  24.0000000 

T3=            10.5000000  10.5000000 

T4=             3.0000000  3.0000000 

Cost  J=       O.9306967E-35  0 . 1177913E-34 

Encounter  Frequency  0.6  Rads  Per  Sec 

Sea  State:            1  2 

Kl=             1.3999996  1.3999996 

Tl=            58.2299957  58.2299957 

T2=            24.0000000  24.0000000 

T3=            10.5000000  10.5000000 

T4=             3.0000000  3.0000000 

Cost  J=       0.1868590E-36  0 . 2364934E-36 

Encounter  Frequency  0.75  Rads  Per  Sec 

Sea  State:            1  2 

Kl=             1.3999996  1.3999996 

Tl=            58.2299957  58.2299957 

T2=            24.0000000  24.0000000 

T3=            10.5000000  10.5000000 

T4=             3.0000000  3.0000000 

Cost  J=       0.1917420E-38  0 . 2426735E-38 

Encounter  Frequency  1.5  Rads  Per  Sec 

Sea  State:            1  2 

Kl=             1.3999996  1.3999996 

Tl=            58.2299957  58.2299957 

T2=            24.0000000  24.0000000 

T3=             10.5000000  10.5000000 

T4=             3.0000000  3.0000000 

Cost  J=       0.1476825E-38  0 . 1869107E-38 
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Sea  State: 

Kl= 

Tl  = 

T2  = 

T3  = 

T4  = 

Cost  J= 

Sea  State: 

Kl= 

Tl= 

T2  = 

T3  = 

T4= 

Cost  J= 

Sea  State: 

Kl= 

Tl  = 

T2  = 

T3  = 

T4  = 

Cost  J= 

Sea  State: 

Kl= 

Tl  = 

T2  = 

T3  = 

T4  = 

Cost  J= 

Sea  State: 

Kl= 

Tl  = 

T2  = 

T3  = 

T4  = 

Cost  J= 

TABLE  5 

Controller  C  For  30°  Encounter  Wave  Angle 

Encounter  Frequency  0.2  Rads  Per  Sec 

1  2 

1,7432022  1.7434111 

35.1575165  35.0142517 

22.0137024  22.0644684 

13.5792084  13.7818375 

22.0728760  22.1337128 

0.594662E+00  0 . 7502608E+00 

Encounter  Frequency  0.4  Rads  Per  Sec 

1  2 

0.2120571  0.1997837 

92.9998932  99.9923096 

24.2121429  24.2589569 

39.9739990  39.9999542 

32.8156433  34.0418396 

0.2301481E-02  0 . 2901770E-02 

Encounter  Frequency  0.6  Rads  Per  Sec 

1.8069019  1.7886524 

1.1657065  1.1832294 

12.4824829  13.3286743 

13.7605896  15.5890350 

9.2525921  9.2333460 

0.1797669E-02  0 . 2275771E-02 

Encounter  Frequency  0.75  Rads  Per  Sec. 

1  2 

2.2712851  2.2712851 

0.84311^a  0.8431141 

24.1657867  24.1657867 

22.2913666  22.2913666 

9.8764191  9.8764191 

0.1094168E-02  0 . 1384705E-02 

Encounter  Frequency  1.5  Rads  Per  Sec 

1  2 

1.6945763  1.6945763 

0.2571034  0.2571034 

25.3036194  25.3036194 

25.4920502  25.4920502 

26.4106903  26.4106903 

0.5423090E-04  0 . 6863636E-04 


24 


Sea  State: 

Kl= 

Tl= 

T2  = 

T3  = 

T4= 

Cost  J= 

Sea  State: 

Kl  = 

Tl= 

T2  = 

T3  = 

T4= 

Cost  J= 

Sea  State: 

Kl= 

Tl= 

T2  = 

T4= 

Cost  J= 

Sea  State: 

Kl= 

Tl= 

T2  = 

T3  = 

T4= 

Cost  J= 

Sea  State: 

Kl  = 

Tl= 

T2  = 

T3  = 

T4  = 

Cost  J= 

TABLE  6 
Controller  C  For  60°  Encounter  Wave  Angle 

Encounter  Frequency  0.2  Rads  Per  Sec 

1  2 

2.1659079  2.1752367 

29.1530609  28.9752655 

18.6070404  18.4342346 

10.7567606  10.5662384 

18.3101044  18.4618378 

0.7797163E=00  0 . 9817770E+00 

Encounter  Frequency  0.4  Rads  Per  Sec 

1  2 

0.8408279  0.8623953 

73.1355896  74.4368439 

53.1565094  57.8289948 

22.6398926  24.3586884 

49.9999390  49.9900360 
0. 78 11032E+000. 978357 lE+00 

Encounter  Frequency  0.6  Rads  Per  Sec 

1  2 

1.0592918  0.9676542 

0.1000010  0.1343303 

19.7355957  25.3748169 

49.9999390  33.4176483 

0.2017351E-02  0 . 2553086E-02 

Ehcounter  Frequency  0.75  Rads  Per  Sec 

1  2 

2.2673655  2.2673655 

0.8431506  0.8431506 

14.1734161  14.1734161 

11.3151855  11.3151855 

15.1259995  15.1259995 

0.2704291E-02  0 . 3422481E-02 

Encounter  Frequency  1.5  Rads  Per  Sec 

1  2 

2.3734770  2.3808193 

0.5810375  0.5764611 

8.2907715  7.0220699 

6.7814789  5.6438484 

7.9672575  9.3380499 

0.1145959E-02  0 . 1450353E-02 
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TABLE  7 

Controller  C  For  90°  Encounter  Wave  Angle 

Encounter  Frequency  0.2  Rads  Per  Sec 

Sea  State:            1  2 

Kl=             1.3473454  1.3487740 

Tl=            35.3837585  35.2085876 

T2=            19.5286560  19.5352631 

T3=            15.7849731  15.8424835 

T4=            19.3971252  19.5299683 

Cost  J=       0.1756000E+01  0 . 2201721E+01 

Encounter  Frequency  0.4  Rads  Per  Sec 

Sea  State:            1  2 

Kl=            1.6298742  1.6293240 

Tl=            31.0220642  31.2911530 

T2=            18.8860321  19.0679169 

T3=            12.3363419  12.2226353 

T4=            18.8753967  19.0280762 

Cost  J=       0.6726980E+00  0 . 8399093E+00 

Encounter  Frequency  0.6  Rads  Per  Sec 

Sea  State:            1  2 

Kl=             1.9138527  1.9098577 

Tl=             1.1895151  1.1946430 

T2=             10.0705585  10.1293569 

T3=             11.7410755  11.8040771 

T4=            10.0737257  10.1994019 

Cost  J=       0.5086016E-01  0 . 6385970E-01 

Encounter  Frequency  0.75  Rads  Per  Sec 

Sea  State*            1  2 

Kl=    '          1.8993120  1.8893242 

Tl=             0.4503812  0.4850104 

T2=            27.3063202  31.5392303 

T3=            21.9325409  27.4768982 

T4=            28.9246674  26.0999603 

Cost  J=        0.2837797E-01  0 . 3586013E-01 

Encounter  Frequency  1.5  Rads  Per  Sec 

Sea  State 

Kl=             2.4287634  2.4257565 

Tl=             0.6342447  0.6373515 

T2=             7.4652443  7.9898472 

T3=              5.9140081  6.1980715 

T4=              5.8591156  5.1957550 

Cost  J=        0.1057035E-01  0 . 1336560E-01 
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TABLE  8 
Controller  C  For  0°  Encounter  Wave  Angle 

Encounter  Frequency  0.2  Rads  Per  Sec 


Sea  State: 
Kl= 
Tl  = 
T2  = 
T3  = 
T4  = 
COST  J=   0, 

4 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

.9972246E-33 

0, 

6 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

,3988896E-32 

0. 

7 

1.3999996 
58.2299957 
24.0000000 
10.5000000 

3.0000000 
,1221600E-31 

Encounter  Fre 

:quency  0.4  Rads 

Per 

Sec 

Sea  State: 
Kl= 
Tl= 
T2= 

III 

COST  J=   0, 

4 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

. 1454213E-34 

0. 

6 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

.5816855E-34 

0. 

7 

1.3999996 
58.2299957 
24.0000000 
10.5000000 

3.0000000 
,1781412E-33 

Encounter  Fre 

iquency  0  .  6  Rads 

Per 

Sec 

Sea  State: 
Kl= 
Tl= 
T2  = 
T3  = 
T4  = 
COST  J=   0 

4 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

.2919672E-36 

0, 

6 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

.1167868E-35 

0. 

7 

1.3999996 
58.2299957 
24.0000000 
10.5000000 

3.0000000 
,3576599E-35 

Encounter  Fre 

:quency  0.75  Rads 

3  Per  Sec 

Sea  State: 
Kl= 
Tl= 
T2  = 
T3  = 
T4= 
COST  J=   0 

4 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

.2995968E-38 

0, 

6 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

.1198388E-37 

0. 

7 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

,3670063E-37 

Encounter  Fre 

tquency  1.5  Rads 

Per 

Sec 

Sea  State: 
Kl= 
Tl  = 
T2  = 
T3  = 
T4= 
COST  J=   0 

4 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

.2307540E-38 

0, 

6 

1.3999996 

58.2299957 

24.0000000 

10.5000000 

3.0000000 

.9230160E-38 

0, 

7 

1.3999996 
58.2299957 
24.0000000 
10.5000000 

3.0000000 
,2826737E-37 
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TABLE  9 
Controller  C  For  30°  Encounter  Wave  Angle 

Encounter  Frequency  0.2  Rads  Per  Sec 


Sea  State 
Kl  = 
Tl  = 
T2  = 
T3  = 
T4  = 
COST  J= 


Sea  State: 
Kl= 
Tl= 
T2  = 
T3  = 
T4= 
COST  J= 


Sea  State 
Kl  = 
Tl  = 
T2  = 
T3  = 
T4  = 
COST  J= 


Sea  State 
Kl  = 
Tl  = 
T2  = 
T3  = 
T4  = 
COST  J= 


1 
35 
22 
13 
22 


4 
7466650 
2978973 
2485657 
6511078 
1172638 


0.9229971E+00 


0 


1,7746067 
36.3333435 
23.2952423 
13.4007721 
23.5081635 
3496342E+01 


1 
48 
45 
25 
23 


8194151 
3066101 
8587494 
3368378 
1087189 


0.9325538E+01 


Encounter  Frequency  0.4  Rads  Per  Sec 


0 
92 
23 
39 
33 


2102537 
8943787 
9096680 
9628601 
1146698 


0 
99 
24 
39 
33 


2002187 
9424896 
3453674 
9452667 
9409027 


0.3595694E-02    0 . 1431070E-01 


0 
95 
24 
38 
33 


2060161 
9424133 
2843170 
9973450 
4644012 


0.4376912E-01 


Encounter  Frequency  0.6  Rads  Per  Sec 


1 

1 

11 

12 

10 


4 
8260603 
1601276 
2977753 
5801849 
2739115 


1 

1 

10 

11 

10 


8175459 
1572495 
3958378 
7078133 
7313347 


0.2808505E-02    0 . 1122463E-01 


1 

1 

10 

11 

10 


8287125 
1652012 
5659571 
7124157 
5683947 


0.3429835E-01 


Encounter  Frequency  0.75  Rads  Per  Sec 


2.2723074 

0.8993683 

22.2714996 

20.4913177 

10.0878286 

0.1709419E-02 


2 

0 

13 

10 

12 


3248615 
8716087 
4890442 
3330994 
5009842 


0.6832119E-02 


2 

0 

13 

10 

13 


3073282 
8750648 
0364380 
1076660 
1554108 


0.2090112E-01 


Encounter  Frequency  1.5  Rads  Per  Sec 


Sea  State 
Kl= 
Tl  = 
T2  = 
T3  = 
T4  = 
COST  J= 


1 

0 

38 

19 

17 


4 
9173651 
6582609 
7910156 
3868103 
9940795 


1 

0 

31 

29 

24 


8726044 
3396787 
8340454 
9999695 
4667816 


0.8474453E-04    0 . 3389677E-03 


1, 

,6678009 

0, 

,3818831 

27, 

,1428680 

26. 

,3686218 

27, 

,3681793 

1038366E-02 
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TABLE  10 
Controller  C  For  60"  Encounter  Wave  Angle 

Encounter  Frequency  0.2  Rads  Per  Sec 


Sea  State 
Kl= 
Tl= 
T2  = 
T3  = 
T4= 
COST  J= 


Sea  State 
Kl= 
Tl= 
T2  = 
T3  = 
T4= 
COST  J= 


Sea  State 
Kl= 
Tl= 
T2  = 
T3  = 
T4= 
COST  J= 


Sea  State 
Kl= 
Tl= 
T2  = 
T3  = 
T4  = 
COST  J= 


2 
29 
18 
10 
18 


4 

1778851 
6119080 
7636871 
2442350 
3449402 


0.1205183E+01 


2.2266541 
31.9857025 
19.4135590 

7.8426695 

19.8131866 

0.4414552E+01 


2 
35 
24 

8 
24 


2981157 
2284546 
2031403 
8281393 
7627716 


0.1080894E+02 


Encounter  Frequency  0.4  Rads  Per  Sec 


0 
77 
57 
21 
58 


4 
8723865 
6901245 
4288635 
6142426 
1057587 


1 

1 

14 

35 

14 


3212109 
1219873 
2327118 
2742767 
2909546 


0.1194492E+01    0 . 3620569E+01 


1.7475338 

0.9073439 

18.7083435 

29.6310730 

20.4921570 

1005195E+02 


Encounter  Frequency  0.6  Rads  Per  Sec 


0 

0 

25 

18 

39 


4 

9356804 
2035089 
2695007 
3522491 
1636200 


0.3151959E-02    0 


1.0065041 

0.1434612 

28.2131348 

16.1927643 

29.4095917 

1260864E-01 


0 

0 

24 

17 

38 


9630518 
1086842 
0041809 
7048645 
1658936 


0.3862068E-01 


Encounter  Frequency  0.75  Rads  Per  Sec 


2 

0 

18 

16 

12 


4 
2310591 
8066239 
4929504 
4596100 
6112747 


2 

0 

14 

12 

14 


2579184 
7918448 
8840103 
1710510 
9709930 


2 

0 

14 

11 

14 


2688522 
7898693 
6062164 
3906097 
7650146 


0.4225172E-02    0 . 1688830E-01 


0.5162057E-01 


Encounter  Frequency  1.5  Rads  Per  Sec 


Sea  State 
Kl  = 
Tl= 
T2  = 
T3  = 
T4= 
COST  J= 


2 

0 
21 
20 

5 


4 
3292980 
5873335 
2806244 
8198700 
3823652 


0.1791609E-02 


2.3723307 
0.5848479 
8.7543087 
7.1684380 
0838490 


0.7157017E-02 


3858051 

6092799 

3835888 

.6424551 

.5094414 


0.2188155E-01 
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TABLE  11 
Controller  C  For  90°  Encounter  Wave  Angle 

Encounter  Frequency  0.2  Rads  Per  Sec 


Sea  State: 
Kl= 
Tl= 
T2  = 
T3  = 
T4= 
COST  J=   0, 

4   ^ 
1.3517475 
35.3476257 
19.6020966 
15.6134491 
19.6564178 
.2690107E+01 

0, 

6 

1.3857212 
36.2312164 
21.4780426 
15.0542984 
21.3916779 
.9237458E+01 

0. 

7 
1.4623213 
42.9548187 
28.4767914 
14.1265106 
28.5155182 
,2028659E+02 

Encounter  Frequency  0.4  Rads 

Per 

Sec 

Sea  State: 
Kl= 
Tl  = 
T2  = 
T3  = 
T4= 
COST  J=   0 

4 
1.0069885 
49.7900391 
36.0238495 
15.7164612 
35.5061340 
.8278778E+01 

0, 

6 

1.1330795 

3.6114454 

18.1521301 

31.1608582 

17.9342651 

. 1712465E+02 

0. 

7 

1.1926146 

0.7214398 

37.0983887 

39.9669037 

38.7577972 

.2298665E+02 

Encounter  Fre 

:quency  0  .  6  Rads 

Per 

Sec 

Sea  state: 
Kl  = 
Tl= 
T2  = 
T3  = 
T4  = 
COST  J=   0 

4 

1.8570919 

1.0337286 

11.9529600 

12.3197098 

11.7377634 

.9655529E-01 

0, 

6 

1.9276266 

0.9287271 

15.0574036 

13.4190674 

15.0851593 

.2932869E+00 

0, 

1.8934374 

0.4949190 

30.3212738 

27.2785339 

26.5594635 

,5224001E+00 

Encounter  Frequency  0.75  Rads  Per  Sec 

Sea  state: 
Kl= 
Tl= 
T2  = 
T3  = 
T4= 
COST  J=   0 

4 

1.9137421 

0.4686716 

26.6723328 

20.3362122 

30.0011749 

.4419509E-01 

0 

6 

1.8589468 

0.3916095 

32.9264069 

27.4372711 

32.4602356 

. 1722615E+00 

0. 

7 

1.7544203 

0.1022463 

42.4169006 

29.9568634 

52.5700073 

,4958326E+00 

Encounter  Frequency  1.5  Rads 

Per 

Sec 

Sea  State: 
Kl= 
Tl  = 
T2  = 
T3  = 
T4  = 
COST  J=   0 

4 
2.3951197 
0.6277218 
4.0210686 
6.0798759 
69.4311066 
.1650761E-01 

0 

6 
2.4414005 
0.6339726 
6.7097816 
4.9043417 
6.7554045 
.6472868E-01 

0. 

7 
2.4639597 
0.6183524 
7.3951139 
5.1939421 
7.3074417 
,  1906425E  +  00 
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V.  IRREGULAR  SEAS -CONTROLLER  DESIGN 

A  sea  state  generator  program  which  genarates  added 
mass,  added  inertia  values  and  in  addition  calculates  the 
forces  and  moments, was  coupled  to  the  fortran  program  shown 
in  Appendix  D,  so  that  the  function  minimization  subroutine 
(BOXPLX)  could  be  used  in  the  presence  of  the  irregular 
sea. The  forces  and  moments  were  stored  in  a  look  up  table 
which  was  coupled  to  the  equation  of  motion. 

The  optimal  gains  obtained  for  0.75  rad/sec  encounter 
frequency  in  the  regular  sea  study  were  used  as  the  initial 
guess  in  order  to  evaluate  the  optimal  controller  parameters 
involving  the  irregular  sea. 

It  is  known  that  sea  is  never  regular  but  actually  is  a 
random  phenomenon  where  waves  are  continually  changing  in 
height,  length  and  breadth. 

Since  the  sea  state  during  this  study  is  represented  by 
irregular  waves  then  the  waves  impinding  on  the  ship  hull 
would  contain  the  total  energy  density  spectrum.  This 
dencity  specrum  is  composed  of  many  frequencies  and  there- 
fore the  response  of  the  ship  would  be  to  an  average  value 
of  added  mass  and  added  inertia. 

For  this  study  where  the  ship's  speed  is  23  knots, the 
controller  has  the  form  of  equation  B.l,  we  used  values  for 
added  mass  and  added  inertia  corresponding  to  0.75  rad/sec 
because  close  to  that  frequency  the  energy  density  is 
maximum. 

The  optimal  controller  parameters  found  are  shown  in 
Table  12 , and  typical  system's  response  is  shown  in  Figures 
5.1,5.2,5.3,5.4,5.5 

These  plots  were  obtained  using  the  program  of  Appendix 
E. 
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A  careful   analysis  of   the  extracted   results  leads   to 
conclude : 

•  The  maximum  deviation  of  controller   parameters  values 
occured  at  30°  encounter  angle. 

•  For  all  encounter   angles  the  maximum  cost   occured  for 
sea  state  7. 

•  For  specific  encounter  angles  the   higher  the  sea  state 
the  higher  the  cost. 
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TABLE  12 
Controller  C  for  Different  Encounter  Wave  Angle 

0°  Encounter  Wave  Angle 
Encounter  Frequency  0.75  Rads  Per  Sec 


Sea  State: 

:       4 

6 

7 

Kl  = 

1.7460003 

1.7460003 

1.7460003 

Tl  = 

35.2969971 

35.2969971 

35.2969971 

T2  = 

22.2480011 

22.2480011 

22.2480011 

T3  = 

13.6510000 

13.6510000 

13.6510000 

T4= 

22.1170044 

22.1170044 

22.1170044 

Cost  J= 

0.9438775E-34 

0.3782326E-34 

0.4605412E-33 

30°  Encounter  Wave  Angle 
Encounter  Frequency  0.75  Rads  Per  Sec 

Sea  State:  4  6  7 

Kl=  2.4596768  1.5688477  2.4804857 

Tl=  88.2797241  47.3260040  56.3383179 

T2=  50.5678864  35.6789203  51.4950714 

T3=  5.2703905  21.5429993  5.7071409 

T4=  95.3189392  25.0237122  91.6153102 

Cost  J=    0.8905333E-01    0 .4461777E-01    0 . 1304857E+00 

60°  Encounter  Wave  Angle 
Encounter  Frequency  0.75  Rads  Per  Sec 

7 

0.2646747 
15.4308896 
35.0908203 
38.4357300 
62.9836578 
0.1193313E-01 

90°  Encounter  Wave  Angle 
Encounter  Frequency  0.75  Rads  Per  Sec 

Sea  State:  4  6  7 

Kl=  0.5902819  2.4912267  2.4547853 

Tl=  87.4059400  55.7961656  32.7171021 

T2=  34.5666444  35.6789999  33.8865433 

T3=  39.6794654  30.2094433  4.9768888 

T4=  25.1668799  25.0789005  25.3325553 
Cost  J=   0.7864684E-01    0 . 1778146E-01    0 . 2522047E+00 


Sea  State: 

4 

6 

Kl  = 

0.3419376 

0.1919854 

Tl  = 

99.6322021 

16.3320312 

T2  = 

35.0456085 

38.0289001 

T3  = 

38.9945677 

34.9145813 

T4= 

25.1149750 

60.8223572 

Cost  J=    0. 

, 1788034E-03 

0.3666982E-02 
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VI.  ADAPTIVE  CONTROL 

A.   STRUCTURE 

When  the  ship  is  moving  in  a  seaway,  the  controller 
parameters  are  changing  due  to  alterations  in  the  sea  state 
and  encounter  angle.  In  addition  we  know  that  using  fixed 
parameters  for  the  controller  over  the  entire  spectrum, it  is 
somewhat  difficult  to  have  an  appropriate  response  of  the 
system. The  adjustment  of  the  controller  parameters  during 
operation  in  a  seaway, can  be  achieved  by  means  of  an  adap- 
tive control. 

The  adaptive  controller  can  be  built  with  digital 
circuits  and  analog  components  as  well.  Analog  system  hard- 
ware has  to  be  designed  for  each  specific  requirement,  and 
any  new  requirement  involves  changes  to  components.  This  is 
a  time  consuming  task.However , f or  simple  control  require- 
ments anolog  systems  still  have  a  possible  economic  advan- 
tage over  digital  systems  [Ref.  11]. 

On  the  other  hand, digital  systems  are  immediately  more 
attractive  when  control  systems  are  required  to  carry  out 
more  and  more  complex  tasks. The  advent  of  microprocessors 
and  associated  components  has  enabled  low  cost  microcom- 
puters to  be  built. These  no  longer  require  special  environ- 
ments and  are  fully  compatible  with  shipboard  use. The 
advantages  of  microcomputer  controls  over  other  types  are 
that  they  are  extremely  reliable  particularly  as  relatively 
fewer  components  are  required  and  hence  they  are  smaller. 
Their  capability  is  greater  than  comparable  analog  systems 
due  to  their  ability  to  carry  out  more  complex  calculations. 
A  major  advantage  is  their  flexibility  while  using  standard 
hardware,    being  able   to  be   reconfigured   for  changes   in 
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system  requirement  without  the  need  to  alter  the  hard- 
ware.This  flexibility  is  achieved  by  the  programmed  software 
which  is  stored  in  the  memory  of  the  computer. 

An  adaptive  control  scheme  is  indicated  in  Figure  6.1 
and  in  Figure  6.2  we  show  analytically  the  components  of 
what  we  call  the  'Decision  Generator' . 

B.   SEARCH  ALGORITHM 

The  parameters  corresponding  to  various  environmental 
conditions (sailing  modes)  are  sorted  in  a  memory. When  the 
sailing  modes  change,  delta(0  )  and  psiCVp)  vary  so  that  the 
value  of  the  performance  index  changes  from  its  theoretical 
value.  This  change  is  detected  by  the  threshold  detector. 

If  the  threshold  is  too  small  the  adaptive  structure 
will  try  new  conditions  when  it  is  not  necessary  and  in  the 
opposite  case  the  adaptive  structure  might  not  adapt  to  new 
sailing  modes. 

Using  the  term  threshold  we  want  to  make  sure  that  only 
significant  changes  in  the  environmental  conditions  will  be 
considered,  otherwise  the  system  will  look  for  a  new  model 
too  often. 

When  the  threshold  detects  that  a  new  set  of  parameters 
should  be  used  the  search  algorithm  will  look  for  conditions 
close  to  the  previous  one  because  the  external  conditions 
are  not  expected  to  vary  rapidly.  Then  a  new  optimal  value 
of  J  is  selected  and  must  be  matched  with  the  performance 
index  computed  for  the  real  psi  and  delta. 

Matching  operation  must  be  done  after  some  delay 
ensuring  that  the  conditions  had  an  effect  on  the  real  cost 
function. 

These  new  parameter  values  will  be  fed  into  the 
controller  and  function  minimization  subroutine  as  well. 
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Figure  6.1    Adaptive  Control  Scheme 

ADVANTAGES 

•  Using  this  adaptive  scheme  we  don't  use  sensors  which 
may  be  unrealistic  and  the  use  of  predetermined  filters 
which  are  so  expensive  can  be  avoided. 
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Figure  6.2    Decision  Generator 

•  Provides  a  good  choice  for  the  function  minimization 
subroutine  to  start  it's  iteration  algorithm  and  there- 
fore we  get  the  optimum  values  as  quickly  as  possible. 

DISADVANTAGES 

•  The  search  algorithm  must  be  carefully  determined. 

•  Since  we  don't  include  sensor's  measurements  there  is 
no   indications  on   how  to   perform   the  search.    This 
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problem  can  be  eliminated  when  the  NAVSTAR/ GLOBAL 
POSITION  SYSTEM  (GPS)  will  be  able  to  provide  precision 
navigation  data. 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS 

A.   CONCLUSIONS 

•  A  properly  designed  controller  can  minimize  the  rudder 
activity  providing  desired  overseas  heading  as  well, and 
therefore  does  result  in  substantial  improvement  in 
propulsion  efficiency. 

•  Actual  savings  in  fuel  cannot  yet  be  determined  since 
there  is  no  information  available  from  the  conventional 
autopilot  and  therefore  there  is  no  possibility  for 
comparison. 

•  It  is  believed  that  the  performance  index  used  in  this 
research  is  a  fairly  adequate  function. Doubts  arise 
from  the  weighting  factor  which  is  included  and  this 
because  lambda  is  based  on  assumptions  and  it's  accu- 
racy is  not  certain. 

•  An  adaptive  controller  which  minimizes  propulsion 
losses  due  to  steering  is  needed  when  environmental 
conditions  and  ship  characteristics  change. 

•  Studying  all  the  investigated  sailing  modes,  it  turns 
out  that  the  cost  surface  is  flat. As  a  consequence 
determination  of  the  controller  parameters  does  not 
require  high  accuracy. 

•  The  study  shows  that  the  use  of  the  third  order  ship 
Nomoto  model  is  a  reasonable  choice  instead  of  using 
the  ship's  equation  of  motion, which  involves  both  the 
sway  and  yaw  equations. 

•  It  can  be  assumed  from  the  work  done  in  this  research 
that  some  of  the  findings  could  have  applicability  to  a 
number  of  ship  types .However  it  is  not  possible  to  make 
inference,   of   a  general   nature  concerning   all  ships 
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from  the  results  of  the  SL-7  which  represents  a  partic- 
ular type  of  a  particular  class  of  ships. 

B.   RECOMMENDATIONS  FOR  FUTURE  STUDY 

•  Additional  work  has  to  be  done  for  obtaining  optimum 
controller  parameters  under  an  expanded  range  of  oper- 
ating conditions.  The  more  optimum  controller  param- 
eter values  available,  the  better  the  determination  of 
the  search  algorithm  recommended  in  the  adaptive 
control  scheme. 

•  A  study  including  the  surge  equation  in  our  ship  model 
is  recommended  since  in  reality  added  resistance  due  to 
steering  reduces  ship's  speed  and  so  far  we  assumed 
constant  speed.  In  addition  with  the  use  of  the  surge 
equation  we  can  determine  actual  energy  losses. 

•  Investigation  of  more  advanced  methods  of  adaptive 
control  based  on  "on-line"  determination  of  plant 
characteristics . 

•  So  far  we  were  interested  in  minimizing  rudder  and 
yawing  activity  to  reduce  propulsion  losses  using  a 
particular  performance  criterion  in  open 
seas . Considering  some  other  types  of  control  such  as 
automatic  control  for  replenishment  requires  an  inves- 
tigation as  to  the  suitability  of  the  cost  function, 
and  a  comparison  with  other  potential  cost  functions. 
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APPENDIX  A 
DETERMINATION  OF  NOMOTO  THIRD  ORDER  MODEL-BOXPLX 

//PROGRA  JOB  (????, 0356) , 'RESEARCH'  ,CLASS  =  G 

//"MAIN  ORG^NPGVMl. ????P 

//  EXEC  F0RTXCG,PARM.F0RT='0PT(2) ' ,IMSL=DP,REGION=1024K 

//FORT. SYS IN  DD  - 

C  THIS  PROGRAM  WILL  OBTAIN  THE  CONTROLLER  OPTIMAL  GAINS. 

C  IT  IS  REFERENCED  IN  CHAPTER  5. 

C  IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS  HAVE  BEEN 

C  OBTAINED  CHANGE  XS(--)  TO  X(''-)  AND  DELETE  XU('-),AND  XL(") 

DIMENSION  XS( 5 ),XU( 5 ),XL(5) 

XS(1)=1. 7466650 

XS(2)=35. 2978973 

XS(3)=22. 2485657 

XS(4)=13. 6511078 

XS(5)=22. 1172638 
C   XS(I)  IS  THE  STARTING  GUESS 

C   XL(I)  IS  THE  LOWER  LIMIT  FOR  THE  I ' TH  VARIABLE 
C   XU(I)  IS  THE  UPPER  LIMIT  FOR  THE  I ' TH  VARIABLE 

XL(1)=0.1 

XU(1)=4.0 

XL(2)=0.1 

XU(2)=80.0 

XL(3)=1.0 

XU(3)=50.0 

XL(4)=1.0 

XU(4)=30.0 

XL(5)=1.0 

XU(5)=50.0 
C   A  DESCRIPTION  OF  THE  FOLLOWING  PARAMETERS 
C   IS  DISCUSSED  IN  BOXPLX 
R=9./13. 
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NTA=1000 
NPR=100 
NAV=0 
NV=5 
IP=0 
C   THE  FOLLOWING  STATEMENT  MUST  BE  CHANGED  TO 
C   CALL  PLANT (X) 
C   IF  ONLY  SIMULATION  IS  WANTED 

CALL  BOXPLX (NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , lER ) 
WRITE  (6,25) 
25      FORMAT (IX,'  OPTIMAL  GAINS ',/ ) 

DO  30  1=1,5 
30      WRITE(6,40)I,XS(I) 
40     F0RMAT(1X, 'X(' ,12, ' )=' ,F14.7) 
STOP 
END 

SUBROUTINE  PLANT (XX) 
C   SUBROUTINE  PLANT (XX)  SIMULATES  THE  SHIP 
COMMON  TDIFF 
REAL" 8  L , L2 , L3 , L4 , L5 , L6 

REAL "8  X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 
REAL" 8  TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 
REAL "8  YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 
REAL" 8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL-8  RHO , IZ , FX , FY ,MZ , XP ,MASS , DELT , MZI , RXI , WA , WE 
REAL" 8  DYAWE , YAWE , YAWC , ISE , ISR , LAMDA , D , RYR , RYI ,MZR 
REAL-8  Kl , Tl , T2 , D , X2 , DX2 , S , RX , RY , RZ , TX , TY , TZ , RXR 
REAL-8  T3,T4,X3,DX3,X4 
DIMENSION  XX(5) 
C 

C   CLOSE  LOOP  ANALYSIS  WITH  FILTER 
C 

C   INITIAL  CONDITIONS  FOR  INTEGRATION 
C   SIMULATION  END  TIME  IN  SECONDS 
ETIME=600.0 
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TIME=0.0 

IC0UNT=1 
C   INITIALIZE  THE  COST  FUNCTION 

ISE=0.0 

ISR=0.0 

TDIFF=0.0 

LAMDA=8.128 
C   GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 

K1=XX(1) 

T1=XX(2) 

T2=XX(3) 

T3=XX(4) 

T4=XX(5) 
C      WRITE(6,1010)  K1,T1,T2 

CIOIO   F0RMAT(1X, 'Kl  =',F15.7,'T1  =',F15.7,'  T2  =',F15.7) 
C   X,XDOT,Y,YDOT  ARE  FIX  COORDINATES  ON  EARTH 

X=0.0 

Y=0.0 

XDOT=0.0 

YDOT=0.0 
C   U,UDOT,V,VDOT  ARE  FIX  COORDINATES  ON  SHIP 

V=0.0 

UDOT=0.0 

VDOT=0.0 

YAW=0.0 

R=0.0 

RDOT=0.0 
C   ORDERED  SPEED  IN   FEET/ SEC 
C    38.82  FT/SEC=23  KNOTS 

UC=38.82 
C   AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

U=UC 
C   D  =  RUDDER  ANGLE 

D=0.0 

L=880.5 
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L2  =  L'-"2 

L3  =  L''L"L 

L4  =  L'*-L3 

L5=L"L4 

L6  =  L''L5 
C   SEA  DISTURBANCE 

C   FORCES  IN  X,Y  DIRECTION  COMPUTED  IN  FORCES 
C   MOMENTS  IN  Z 

FX=0. 

FY=0. 

MZ  =  0. 

RXR=-0.15744D+05 

RXI=-0.19950D+06 

RYR=0.52365D+04 

RYI=0.18699D+06 

MZR=-0.29870D+08 

MZI=-0.35751D+07 
C     RXR=-0.50230D+04 
C      RXI=0.12712D+05 
C     RYR=0.35290D+04 
C     RYI=-0.31909D+05 
C     MZR=0.38826D+07 
C      MZI=-0.64313D+07 
C      RXR=0.28540D+04 
C     RXI=-0.99574D+04 
C     RYR=-0.85441D+04 
C      RYI=0.39595D+05 
C     MZR=-0. 13014D+08 
C     MZI=0.11348D+08 
C      RXR=-0.75642D+04 
C     RXI=0.83497D+04 
C     RYR=0.23379D+05 
C     RYI=-0.81502D+05 
C     MZR=0.28622D+07 
C     MZI=-0.19388D+08 
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C  RXR=-0.37916D+04 

C  RXI=0.16381D+04 

C  RYR=-0.76647D+05 

C  RYI=-0.37685D+05 

C  MZR=-0.83915D+07 

C  MZI=-0.53176D+07 

RX  =  D S QRT  ( RXR '- "  2  +  RXI  -' -  2  ) 

RY=DSQRT  (RYR'"''2  +  RYI--2 ) 

RZ^DSQRT  (MZR-'"''2  +  MZI" "2 ) 

TX=DATAN(RXI/RXR) 

TY=DATAN(RYI/RYR) 

TZ=DATAN(MZI/MZR) 
C   SIGNIFICANT  WAVE  HEIGHT;  SEA  STATE  1-5,2-10,3-15 
C   4-17.5,5-22.5,6-27,7-35,8-42,9-60 

WA=17.5 
C   ENCOUNTER  FREQUENCY  .1,. 2, .3, .4, .5, .6, .75, 1.0, 1.5, 2. 5 

WE=0.2 
C   HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  HERE  AS 
c   PARAMETERS 

RHO=1.9876 

MASS  =  (  .  0044 )  "  (  .  5 ''RH0''-L3  ) 

IZ= (0 . 00028 )"( . 5-RHO-L5 ) 

YAWE=0.0 

X2=0.0 

DX2=0.0 

X3=0.0 

DX3=0.0 
200  CONTINUE 

S  =  DSQRT  (u-*- '''2+V'-"2  ) 
C   INPUT  YAW  COMMAND 

YAWC=0.0 

IF  (TIME. GE. 0.0)  YAWC=0.0 
C   ERROR  SIGNAL  TO  DRIVE  RUDDER (YAW  ACTUAL  -  YAW  ORDERED) 
C   (  COMPENSATOR  FILTER  ) 

YAWE=YAW  -  YAWC 
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DX2=(YAWE-X2)/T2 

X4=K1"(T1"DX2+X2) 

DX3=(X4-X3)/T4 

D=(T3"DX3+X4) 
C   AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 
C   XUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED 
C   FOR  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  , 
C   ENCOUNTER  FREQUENCY 

XUDOT = ( - . 0001 ) " ( • 5 "RH0-L3 ) 

XU=(-0.0253)^M.5"RHO''-L2"S) 

XUU=  (  -  0  .  0003  )  '•- ( .  5  "RH0^'^L2 ) 

XVR= ( 0 . 00  3  9 ) " ( . 5 "RHO "L3 ) 

XVV= ( - . 00 12 ) " ( . 5 "RH0-L2 ) 

XDD=  ( -  0  .  0005  ) "  (  .  5  "RH0-L2 "S''"2 ) 
C   LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 
C     YV=  ( - 0  .  00758  ) "  ( .  5^"RH0^'^L2^*^S ) 

YR= (0 . 0023 )"( . 5-RH0-L3-S) 

YD=  (0  .  00145  )"(  .  5-RH0-L2-S  — 2  ) 

YVVR= ( 0 . 0 1 ) " ( . 5 " RHO "L3 / S ) 

YVRR= ( -  0 . 0  0  8 ) " ( . 5 -RHO - L4 / S ) 

YVVV= ( - 0 . 0  3 ) " ( . 5 "RHO "L2 / S ) 

YRRR= ( 0 . 0  0  3 ) " ( . 5  - RHO - L5 / S ) 

YDDD=  (  -  0  .  0005  )  "  (  •  5  "RH0-L2 - S - -'2  ) 
C   YUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE 
C   CHANGED  FOR  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  , 
C   ENCOUNTER  FREQUENCY 
C      YVDOT= ( -  0 . 003  9 ) " ( . 5 "RH0-L3 ) 

C   SPEED=23  KNOTS,  ENCOUNTER  ANGLE  =  30  ,  ENCOUNTER 
C   FREQUENCY  =0.2 

YVDOT=-0.30908D+07 

YV=-0.81271D+04 
C      YVDOT=-0.36185D+07 
C      YV=-0.24757D+06 
C      YVDOT=-0.32890D+07 
C      YV=-0.11775D+07 
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C     YVDOT=-0.23038D+07 

C     YV=-0.18267D+07 

C     YVDOT=-0.59800D+06 

C     YV=-0.13260D+07 

C   MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 

NV=  ( -0  .  00213  )"  ( .  5^"RHO^-L3"S ) 
C     NR=  (  -  0  .  00 105  )  "  (  .  5  "RH0^'^L4"S ) 

ND=  (  -0  .  0007  )"(  .  5"RH0"L3'''S""2 ) 

NVVR=  ( -  0  .  0 15  )  "  ( .  5  '-RHO "LA /  S ) 

NVRR=  (  -  0  .  0 0 8  )  "  (  .  5  -RHO --L5  /  S  ) 

NVVV=  (  0  .  0 1 )  "  (  .  5  "RHO "L3  /  S  ) 

NRRR=  ( -  0  .  006  )  '•-  (  .  5  "RH0'-L6  /  S ) 

NDDD=  ( 0  .  000 1 )  '•- (  •  5 "RH0'''L3 " S ""2  ) 
C   NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE 
C   CHANGED  FOR  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  , 
C   ENCOUNTER  FREQUENCY 
C 

C      NRDOT=  (  -  0  .  00027  ) -'  (  .  5 --RH0-L5  ) 

C   SPEED=23  KNOTS,  ENCOUNTER  ANGLE  =  90  ,  ENCOUNTER 
C   FREQUENCY  =0.2 

NRDOT=-0.26251D+12 

NR=-0.53637D+09 
C  NRDOT=-0.20125D+12 

C  NR=-0.94970D+10 

C  NRDOT=-0.18671D+12 

C  NR=-0.46860D+11 

C  NRDOT=-0.14518D+12 

C  NR=-0.87538D+11 

C  NRDOT=-0.37261D+11 

C  NR=-0.69856D+11 

C      REGULAR  WAVE    SEA   STATE 

FX=WA"RX"DCOS(WE"TIME+TX) 

FY=WA"RY"DCOS (WE-TIME+TY ) 

MZ=WA"RZ"DCOS (WE-TIME+TZ ) 
C   U  ACTUAL  SPEED 
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C   UC  COMMANDED  SPEED 
C   XP  =  PROPELLER  THRUST 

XP=-XUU"UC""2 
C   EQUATIONS  OF  MOTION 

C     UDOT=(  (XVR  +  MASS)"V"R  +  XUU''-U-''"2  +  XW-V  —  Z 
C     1  +  XDD''-D"D  +   FX  +  XP  )  /  (MASS-XUDOT) 

VDOT=(YV"V  +  (YR-MASS''-U)"R  +  YD-D  +  YVVR-V  —  Z-R 

1  +YVRR"V"R"'''2  +YRRR"R"'-3 

2  +  YDDD'-D''-"3  +  FY  )  /  (MASS-YVDOT) 
RDOT=(  NV-V  +  NR-R  +  ND-D  +  NVVR"V""2"R 

1  +  NVRR"V"R""2  +NVVV"V""3 

2  +  NRRR''-R""3  +  NDDD-D—O  +  MZ  )/(IZ-NRDOT) 
C   WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.il)  GO  TO  50 
GO  TO  300 
C   CONVERT  RADIANS  TO  DEGREES 
50   YAWDEG=  YAW -5 7. 29 6 
RDEG=R"5  7.296 
RDDEG  =  RDOT'-57.296 
DDEG=D"57.296 
YAWC=YAWC"57.296 
C      WRITE  (6,100)  TIME,XP,X,XDOT,Y,YDOT 
C     1   ,UC,U,UDOT,V,VDOT,YAWC,YAWDEG,RDEG,RDDEG,DDEG 
100  F0RMAT(1X, 'TIME=' ,F8.3, '  SEC   XP=',F10.2,'  LBF 
1  X=' ,F8.2, 

1'  FT   XDOT=' ,F8.4, '  FT/ SEC   Y=',F8.2,'  FT   YDOT= ' 
1  ,F8.4,'  FT/SEC 

l',/,2X,'  UC=',F8.4,'  FT/SEC    U=',F8.4,'  FT/SEC 
1  UDOT=' ,F10.6, 

1  '  FT/ SEC" "2    V=',F8.4,'  FT/ SEC    VDOT= ' , FIO . 6 , 
1  '  FT/ SEC" -2' 

1  ,/,2X, 'YAWC=' ,F8.4, '  DEG   YAW= ' 
1    ,F15.7,'  DEG    YAW  RATE=' ,F15.7, '  DEG/ SEC 
1  YAW  ACCEL=' ,F15.7, ' 
1  DEG/ SEC" -2' ,/ ,2X, 'RUDDER  =',F15.7/  DEG  ',/) 
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IC0UNT=1 
C   TEST  IF  WANT  TO  STOP 

300   IF  (TIME.GE.ETIME)  GO  TO  400 
C   INTEGRATION  STEP  SIZE  CELT 

DELT^l.O 
C   INTEGRATION 

U=U+UDOT'*-DELT 

V=V+VDOT"DELT 

R=R+RDOT--DELT 

YAW=YAW+R"DELT 

X2=X2+DX2"DELT 

X3=X3+DX3"DELT 
C   CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 
C      XDOT=U"DCOS(YAW)-V"DSIN(YAW) 
C      YDOT=U"DSIN(YAW)+V"DCOS(YAW) 
C      X=X+XDOT"DELT 
C      Y=Y+YDOT"DELT 

TIME=TIME+DELT 

IC0UNT=IC0UNT+1 

ISE  =  ISE  +  LAMDA'-YAWE"'-2 

ISR=ISR  +  D""2 

GO  TO  200 
C   J=TDIFF=  COST  FUNCTION 
400   TDIFF=ISE+ISR 

WRITE (6, 500)  ISE,ISR,TDIFF,K1,T1,T2,T3,T4 
500   FORMAT ('  ' , IX, ' TOTAL= ' , F15 . 7 , 2X, 

1  'Kl=' ,F15.7,2X, 'Tl=' ,F15.7,2X, 'T2=' ,F15.7,2X, 
1  'T3=' ,F15.7,2X, 'T4=' ,F15.7) 

RETURN 

END 

SUBROUTINE  BOXPLX  (NV,NAV ,NPR,NTZ ,RZ ,XS , IP ,BU ,BL, 
1YMN,IER) 
C 

DIMENSION  V(50, 50),  FUN(50),  SUM(25),  CEN(25), 
IXS (NV ) , BU (NV ) , BL (NV ) 
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c 
c 


c 
c 
c 
c 


KV  =  5 

EP  =  l.E-6 

NTA  =  2000 

IF  (NTZ.GT.O)  NTA  =  NTZ 

R  =  RZ 

IF  (R.LE.O. .OR.R.GE.l. )  R=l./3. 

NVT  =  NV+NAV 

TOTAL  VARS,  EXPLICIT  PLUS  IMPLICIT 
NT  =  0 

CURRENT  TRIAL  NO. 
NPT  =  0 

CURRENT  NO.  OF  PERMISSIBLE  TRIALS 
NTFS  =  0 
CURRENT  NO.  OF  TIMES  F  HAS  BEEN  ALMOST  UNCHANGED 

CHECK  FEASIBILITY  OF  START  POINT 

DO  4  1=1, NV 

VT  =  XS(I) 

IF  (BL(I).LE.VT)  GO  TO  1 

II  =  -I 

VT  =  BL(I) 

GO  TO  2 

1  IF  (BU(I).GE.VT)  GO  TO  3 
II  =  I 

VT  =  BU(I) 

2  IF  (NPR.GT.O)  WRITE  (6,49)  II 

3  V(I,1)  =  VT 
CEN(I)  =  VT 

IF  (IP.EQ.l)  GO  TO  4 

BL(I)  =  BL(I)+AMAX1(EP,EP"ABS(BL(I))) 

BU(I)  =  BU(I)-AMAX1(EP,EP"ABS(BU(I))) 

4  SUM(I)  =  VT 
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c 
c 

NCE  =  1 
C  NUMBER  OF  CONSTRAINT  EVALUATIONS 

1=1 

IF  (KE(V(l,l)).EQ.O)  GO  TO  5 
IF  (NPR.LE.O)  GO  TO  12 
WRITE  (6,50) 
GO  TO  12 

5  NFE  =  1 
C 

C    NUMBER  OF  VERTICES  (K)  =  2  TIMES  NO.  OF  VARIABLES 

K  =  2-NV 
C 
C    NUMBER  OF  DISPLACEMENTS  ALLOWED. 

NLIM  =  5-'-NV+10 
C 

C    NUMBER  OF  CONSECUTIVE  TRIALS  WITH  UNCHANGED 
c    FE  TO  TERMINATE. 

NCT  =  NLIM+NV 

ALPHA  =  1.3 

FK  =  K 

FKM  =  FK-1. 

BETA  =  ALPHA+1. 
C 
C    INSURE  SEED  OF  RANDOM  NUMBER  GENERATOR  IS  ODD. 

IQR  =  R-1.E7 

IF  (MOD(IQR,2).EQ.O)  IQR=IQR+101 
C 
C  SET  UP  INITIAL  VERTICES 

FUN(l)  =  FE(V(1,1)) 

YMN  =  FUN(l) 

6  FI  =  1. 
FUNOLD  =  FUN(l) 
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DO  15  1=2, K 

FI  =    FI+1. 

LIMT  =  0 

7  LIMT  =  LIMT+1 

c 

c 

END  CALCULATION 

IF  (LIMT.GE.NLIM)  GO  TO  11 
C 

DO  8  J=1,NV 
C 
C    RANDOM  NUMBER  GENERATOR   (RANDU) 

IQR  =  IQR-65539 

IF  (IQR.LT.O)  IQR  =  IQR+2147483647+ 1 

RQX  =  IQR 

RQX    =    RQX-.4656613E-9 

V(J,I)  =  BL(J)+RQX"(BU(J)-BL(J)) 

IF  (IP.EQ.l)  V(J,I)=AINT(V(J,I)+.5) 

8  CONTINUE 
C 

DO  10  L=1,NLIM 
NCE  =  NCE+1 

IF  (KE(V(l,I)).EQ.O)  GO  TO  13 
C 

DO  9  J=1,NV 

VT  =  .5"(V(J,I)+CEN(J)) 

IF  (IP.EQ.l)  VT  =  AINT(VT+.5) 

V(J,I)  =  VT 

9  CONTINUE 
C 

10  CONTINUE 
C 

11  IF  (NPR.LE.O)  GO  TO  12 
WRITE  (6,51)  I 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,I,FUN,CEN,I) 

12  lER  =  -1 
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GO  TO  48 
C 

13  DO  14  J=1,NV 

SUM(J)  =  SUM(J)+V(J,I) 

14  CEN(J)  =  SUM(J)/FI 
C 

C    TRY  TO  ASSURE  FEASIBLE  CENTROID  FOR  STARTING. 
NCE  =  NCE+1 

IF  (KE(CEN).EQ.O)  GO  TO  60 
SUM(J)  =  SUM(J)  -V(J,I) 
GO  TO  7 
60  NFE  =  NFE+1 

FUN(I)  =  FE(V(1,I)) 

15  CONTINUE 
C 

C    END  OF  LOOP  SETTING  OF  INITIAL  COMPLEX. 

IF  (NPR.LE.O)  GO  TO  17 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,0) 
C 
C    FIND  THE-  WORST  VERTEX,  THE  'J'TH. 

J  =  1 
C 

DO  16  1=2, K 

IF  (FUN(J) .GE.FUN(I))  GO  TO  16 

J  =  I 

16  CONTINUE 
C 

C    BASIC  LOOP.   ELIMINATE  EACH  WORST  VERTEX 
C    IN  TURN.  IT  MUST  BECOME  NO  LONGER  WORST, NOT 
C    MERELY  IMPROVED.   FIND  NEXT-TO-WORST  VERTEX, 
C    THE  'JN'TH  ONE. 

17  JN  =  1 

IF  (J.EQ.l)  JN  =  2 
C 

DO  18  1=1, K 
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IF  (I.EQ.J)  GO  TO  18 

IF  (FUN(JN) .GE.FUN(I))  GO  TO  18 

JN  =  I 

18  CONTINUE 
C 

C    LIMT  =  NUMBER  OF  MOVES  DURING  THIS  TRIAL  TOWARD 
C     THE  CENTRIOD  DUE  TO  FUNCTION  VALUE. 

LIMT  =  1 
C 

C    COMPUTE  CENTROID  AND  OVER  REFLECT  WORST  VERTEX. 
C 

DO  19  1=1, NV 

VT  =  V(I,J) 

SUM(I)  =  SUM(I)-VT 

CEN(I)  =  SUM(I)/FKM 

VT  =  BETA"CEN(I)-ALPHA"VT 

IF  (IP.EQ.l)  VT  =  AINT(VT+.5) 
C 
C    INSURE  THE  EXPLICIT  CONSTRAINTS  ARE  OBSERVED. 

19  V(I,J)  =  AMAX1(AMIN1(VT,BU(I)),BL(I)) 
C 

NT  =  NT+1 
C 

C    CHECK  FOR  IMPLICIT  CONSTRAINT  VIOLATION. 
C 

20  DO  25  N=1,NLIM 
NCE  =  NCE+1 

IF  (KE(V(1,J)) .EQ.O)  GO  TO  26 
C 

C    EVERY  'KV'TH  TIME,  OVER-REFLECT  THE  OFFENDING 
C    VERTEX  THROUGH  THE  BEST  VERTEX. 

IF  (MOD(N,KV) ,NE.O)  GO  TO  22 

CALL  FBV  (K,FUN,M) 
C 

DO  21  1=1, NV 
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VT  =    BETA"V(I,M)-ALPHA"V(I,J) 
IF  (IP.EQ.l)  VT  =  AINT(VT+.5) 

21  V(I,J)  =    AMAX1(AMIN1(VT,BU(I)) ,BL(I)) 
C 

GO  TO  24 
C 

C    CONSTRAINT  VIOLATION:   MOVE  NEW  POINT  TOWARD  CENTROID, 
C 

22  DO  23  1=1, NV 

VT  =  .5"(CEN(I)+V(I,J)) 

IF  (IP.EQ.l)  VT  =  AINT(VT+.5) 

V(I,J)  =  VT 

23  CONTINUE 
C 

24  NT  =    NT+1 

25  CONTINUE 
C 

lER  =  1 
C 

C    CANNOT  GET  FEASIBLE  VERTEX  BY  MOVING  TOWARD  CENTROID, 
C    OR  BY  OVER-REFLECTING  THRU  THE  BEST  VERTEX. 

IF  (NPR.LE.O)  GO  TO  42 

WRITE  (6,52)  NT, J 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,J) 

GO  TO  42 
C 
C    FEASIBLE  VERTEX  FOUND , EVALUATE  THE  OBJECTIVE  FUNCTION 

26  NFE  =  NFE+1 
FUNTRY  =  FE(V(1,J)) 

C 

C    TEST  TO  SEE  IF  FUNCTION  VALUE  HAS  NOT  CHANGED. 

AFO  =  ABS( FUNTRY -FUNOLD) 

AMX  =  AMAX 1(ABS(EP "FUNOLD ) ,EP) 
C 
C    ACTIVATE  THE  FOLLOWING  TWO  STATEMENTS 
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C    FOR  DIAGNOSTICS  PURPOSES  ONLY. 

C      WRITE  (6,99)  J ,AFO ,AMX,FUNTRY,FUNOLD,FUN( J) , 

C     IFUN(JN) ,NTFS,N 

C   99  FORMAT  ( IX , 13 , 6E15 . 7 , 215 ) 

IF  (AFO.GT.AMX)  GO  TO  27 

NTFS  =  NTFS+1 

IF  (NTFS.LT.NCT)  GO  TO  28 

lER  =  0 

IF  (NPR.LE.O)  GO  TO  42 

WRITE  (6,53)  K 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,0) 

GO  TO  42 

27  NTFS  =  0 
C 

C    IS  THE  NEW  VERTEX  NO  LONGER  WORST? 

28  IF  (FUNTRY.LT.FUN(JN))  GO  TO  34 
C 

C    TRIAL  VERTEX  IS  STILL  WORST;  ADJUST  TOWARD  CENTROID 
C    EVERY  'KV'TH  TIME,  OVER-REFLECT  THE  OFFENDING 
C    VERTEX  THROUGH  THE  BEST  VERTEX. 

LIMT  =  LIMT+1 

IF  (MOD (LIMT, KV ) .NE.O)  GO  TO  30 

CALL  FBV  (K,FUN,M) 
C 

DO  29  1=1, NV 

VT  =  BETA"V(I,M)-ALPHA'-V(I,J) 

IF  (IP.EQ.l)  VT  =  AINT(VT+.5) 

29  V(I,J)  =  AMAX1(AMIN1(VT,BU(I)) ,BL(I)) 
C 

GO  TO  32 
C 

30  DO  31  1=1, NV 

VT  =  .5"(CEN(I)+V(I,J)) 

IF  (IP.EQ.l)  VT  =  AINT(VT+.5) 

V(I,J)  =  VT 
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31  CONTINUE 
C 

32  IF  (LIMT.LT.NLIM)  GO  TO  33 
C 

C    CANNOT  MAKE  THE  ' J ' TH  VERTEX  NO  LONGER  WORST 
C    BY  DISPLACING  TOWARD  OVER-REFLECTING 
C    THRU  THE  BEST  VERTEX. 

lER  =  2 

IF  (NPR  .LE.  0)   GO  TO  42 

WRITE  (6,52)   NT,  J 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN, J) 

GO  TO  42 

33  NT  =  NT+1 
GO  TO  20 

C 

C    SUCCESS:   WE  HAVE  A  REPLACEMENT  FOR  VERTEX  J. 

34  FUN(J)  =  FUNTRY 
FUNOLD  =  FUNTRY 
NPT  =  NPT+1 

C 

C    EVERY  100 'TH  PERMISSIBLE  TRIAL,  RECOMPUTE 

C    CENTRIOD  SUMMATION  TO  AVOID  CREEPING  ERROR. 

IF  (MOD(NPT,100) .NE.O)  GO  TO  37 
C 

DO  36  1=1, NV 

SUM(I)  =  0. 
C 

DO  35  N=1,K 

35  SUM(I)  =  SUM(I)+V(I,N) 
C 

CEN(I)  =  SUM(I)/FK 

36  CONTINUE 
C 

LC  =  0 
GO  TO  39 
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c 
c 


37  DO  38  1=1, NV 

38  SUM(I)  =  SUM(I)+V(I,J) 

LC  =  J 

39  IF  (NPR.LE.O)  GO  TO  40 

IF  (MOD(NPT,NPR) .NE.O)  GO  TO  40 


CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,LC) 
C 

C   HAS  THE  MAX.  NUMBER  OF  TRIALS  BEEN  REACHED 
C    WITHOUT  CONVERGENCE? 
C    IF  NOT,  GO  TO  NEW  TRIAL. 

40  IF  (NT.GE.NTA)  GO  TO  41 
C 

C    NEXT- TO- WORST  VERTEX  NOW  BECOMES  WORST. 
J  =  JN 
GO  TO  17 

41  lER  =  3 

IF  (NPR.GT.O)  WRITE  (6,54) 
C 

C    COLLECTOR  POINT  FOR  ALL  ENDINGS. 

C  1)   CANNOT  DEVELOP  FEASIBLE  VERTEX.  lER  =  1 

C  2)  CANNOT  DEVELOP  A  NO-LONGER-WORST  VERTEX. lER  =  2 
C  3)  FUNCTION  VALUE  UNCHANGED  FOR  K  TRIALS.  lER  =  0 
C  4)   LIMIT  ON  TRIALS  REACHED.  lER  =  3 

C  5)   CANNOT  FIND  FEASIBLE  VERTEX  AT  START.    lER  =  -1 

42  CONTINUE 
C 

C    FIND  BEST  VERTEX. 

CALL  FBV  (K,FUN,M) 

IF  (IER.GE.3)  GO  TO  44 
C 
C    RESTART  IF  THIS  SOLUTION  IS  SIGNIFICANTLY  BETTER 
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C    THAN  THE  PREVIOUS,  OR  IF  THIS  IS  THE  FIRST  TRY. 
IF  (NPR.LE.O)  GO  TO  43 
WRITE  (6,55)  (M,YMN,FUN(M)) 

43  IF  (FUN(M) .GE.YMN)  GO  TO  47 

IF  (ABS(FUN(M)-YMN).LE.AMAX1(EP,EP'-YMN))  GO  TO  47 
C 
C    GIVE  IT  ANOTHER  TRY  UNLESS  LIMIT  ON  TRIALS  REACHED. 

44  YMN  =  FUN(M) 
FUN(l)  =  FUN(M) 

C 

DO  45  1=1, NV 
CEN(I)  =  V(I,M) 
SUM(I)  =  V(I,M) 

45  V(I,1)  =   V(I,M) 
C 

DO  46  1=1, NVT 

46  XS(I)  =  V(I,M) 
C 

IF  (IER.LT.3)  GO  TO  6 

47  IF  (NPR.LE.O)  GO  TO  48 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,V(1,M) ,-1) 
WRITE  (6,56)  FUN(M) 

48  RETURN 
C 

49  FORMAT  (50H0INDEX  AND  DIRECTION  OF 
lOUTLYING  VARIABLE  AT  STARTI5) 

50  FORMAT  (50H0IMPLICIT  CONSTRAINT 
IVIOLATED  AT  START.  DEAD  END.) 

51  FORMAT  ('OCANNOT  FIND  FEASIBLE ', 14 ,' TH 
IVERTEX  OR  CENTROID  AT  START.') 

52  FORMAT  ( lOHOAT  TRIAL  I4,54H  CANNOT  FIND 
IFEASIBLE  VERTEX  WHICH  IS  NO  LONGER 
1W0RST,I4,15X, 'RESTART  FROM  BEST  VERTEX.') 

53  FORMAT  (40H0FUNCTION  HAS  BEEN  ALMOST 
lUNCHANGED  FOR  15, 7H  TRIALS) 
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54  FORMAT  (27H0LIMIT  ON  TRIALS  EXCEEDED.  ) 

55  FORMAT  (' OBEST  VERTEX  IS  NO. ',13,' 

1  OLD  MIN  WAS  'E15.7',   NEW  MIN  IS  ',E15.7) 

56  FORMAT  ( ' OMIN  OBJECTIVE  FUNCTION  IS  ',E15.7) 
END 

SUBROUTINE  FBV  (K,FUN,M) 
DIMENSION  FUN (50) 
M  =  1 

DO  1  1=2, K 

IF  (FUN(M).LE.FUN(I))  GO  TO  1 
M  =  I 
1  CONTINUE 

RETURN 

END 

SUBROUTINE  BOUT  (NT ,NPT ,NFE ,NCE ,NV,NVT ,V ,K ,FN , C , IK) 

DIMENSION  V(50, 50) ,  FN(50),  C(25) 

WRITE  (6,4)  NT,NPT,NFE,NCE 

DO  1  1=1, K 

WRITE  (6,5)  FN(I), (V(J,I) ,J=1,NV) 

IF  (NVT.LE.NV)  GO  TO  1 

NVP  =  NV+1 

WRITE  (6,6)  (V(J,I),J=NVP,NVT) 

1  CONTINUE 

IF  (IK.NE.O)  GO  TO  2 

WRITE  (6,7)  (C(I),I=1,NV) 
RETURN 

2  IF  (IK.GE.O)  GO  TO  3 
WRITE  (6,8)  (C(I),I=1,NV) 
RETURN 

3  WRITE  (6,9)  IK, (C(I) ,I=1,NV) 
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5 


RETURN 
C 

4  FORMAT  ('ONO.  TOTAL  TRIALS  =  ',I5,4X 
I'NO.  FEASIBLE  TRIALS  =    ',I5,4X,'NO.  FUNCTION 
lEVALUATIONS  =  ',I5,4X,'N0.  CONSTRAINT  EVALUATI 
IONS  =  ',I5/'0      FUNCTION  VALUE' ,6X, 'INDEPENDENT 
IVARIABLES/DEPENDENT  OR  IMPLICIT  CONSTRAINTS') 

5  FORMAT  (IH  , E18 . 7 , 2X , 7E14 . 7/ (21X , 7E14 . 7 ) ) 

6  FORMAT  (2IX,7E14.7) 

7  FORMAT  (lOHOCENTROID  IIX, 7E14 . 7/ (21X, 7E14 . 7 ) ) 

8  FORMAT  ('0   BEST  VERTEX' , 7X, 7E14 . 7/ (21X , 7E14 . 7 ) ) 

9  FORMAT  ('OCENTROID  LESS  VX' , 12 , 2X, 7E14 . 7/ 
1(21X,7E14.7)) 

END 

FUNCTION  FE(X) 

DIMENSION  X(5) 

COMMON  TDIFF 

CALL  PLANT (X) 

FE=TDIFF 

RETURN 

END 

FUNCTION  KE(X) 

DIMENSION  X(5) 

KE  =  0 

RETURN 

END 
//GO.SYSIN  DD  '- 
/* 
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APPENDIX  B 
CONTROLLER  DESIGN 


A.   CONTROLLER  "C 


Al 


To 


I  _ 

I 

I 

IDX' 


1 

T9 


►  K^Tj 


An 


1 

T, 


T, 


*^ 


DX' 


Figure  B.l    Block  Diagram  of  Controller  C 
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Figure  B.l   corresponds  to  controller   'C'  which  has  the 


form 


K^d  +  T^S)  (l+T^S) 
(1+T23)  (I+T3S) 


(B.l) 


Verifying  that  equation  B.l   corresponds  to  controller 
'C'  we  have: 


A— -— -+K1T1DX2 

^  l4TnS    ^      ^ 


(B.2) 


DX2   s 


(B.3) 


X2   1 


X9  = 


X 


1  +  ToS 


X,-1+T2S     ^  I+V2 
1  (B.4) 

Substituting  equations  B.3  ,   B.4  into  equation  B.2  we 

have : 


XiKi     K2TiX-,S 


I+T2S     I+T2S 
We  also  have: 


(B.5) 


Xu= — +  T^DX3 

"^"^    I+T3S     ^   -^ 


(B.6) 


DX3"  S 


DX3=X3S 


(B.7) 
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X 


3 L_     X  —  A  (B.8) 


Substituting  equations  B.7  ,   B.8  into  equation  B.6  we 
have : 

A    ^  T^AS 
^^"1+138    I+T3S 

(B.9) 
Now   substituting  equation   B.5  into   equation  B.9   we 


have 


XiKi  XiKiTiS 


(1+T2S)(1+T3S)     (l+T2S)(l+ToS) 

X^J^lTi^S  X^KjT^TijS 

(1+T23)(1+T3S)     (1+T2S)(1+T3S)  (B.IO) 


Finally   rearranging    terms   equation   B.IO    can   be 
written: 


X^_Kl(l^TlS-^T„S-^T;^T^S^)  _  K^  (l^T^S  )  (1^T^S  ) 

Xi    (1+T2S)(1+T3S)         C1+T2S)(1+T3S)        (B.ll) 


B.   CONTROLLER  "B" 

Figure  B.2  corresponds   to  controller  'B'  which   has  the 
form: 
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X. 


T' 


-^t 


DX. 


I 

To 


4  DX  -, 


T- 


X 


Figure  B.2    Block  Diagram  of  Controller  B 


Ki(l+TiS) 


(1+T2S)(1+T3S) 


(B.12) 

Verifying  that  equation  B.12  corresponds  to  Controller 
' B '  we  have : 


A  = 


KtX 


I'^l 


iTtTs  "^  ^lTiDX2 


(B.13) 
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X2   _J_ 

DXT'  S 


DX2=  X2S  (B.14) 


X2   1  Xi 

X 


-=— ^ X2  =  i.T  Q  (B.15) 

3^  T2S+I       l  +  T^S 


Substituting  equations  B . 14  ,   B.15  into  equation  B.13 
we  have : 


11  ^"jT-iX-iS 


A  = +. 


I+T2S  l+T^S  (B.16) 

In  a  similar  way  we  also  can  derive  equation  B.17 


X3=    ' 


I  +  T3S 


(B.17) 
Substituting  equation  B.16  into  equation  B.17  we  have 


KnXn  K-iTiXtS 

X3=  i— i +       -L   -t  -L 


(1+T2S)(1+T3S)    (1+T2S)(1+T3S)  (B.18) 

Finally  rearranging  terms  equation  B.18  becomes: 


X3   K^d  +  TiS) 

x7^(l+ToS)(l+ToS)  (B.19) 
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C.   CONTROLLER  "A' 


p 

►   K^T^      

• 

y-l  ^ 

1 

^2 

-Oi' 

1 

Ml         __ 

C 

X2                            +      j^^  X3 

'DXj 

1 

T2     ^ 

h-J 

Figure  B.3    Block  Diagram  of  Controller  A 

Figure  B.3  corresponds   to  controller  'A'  which   has  the 
form: 


K^Cl+T^S) 


I  +  T2S 


(B.20) 
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Verifying  that  equation  B.20  corresponds  to  Controller 
'A'  we  have: 


X,K, 


X3=— +  K1T1DX2 

^    1+ToS     -L  J-   ^ 


But  since  we  know  that 


DXo   1 

1= DX2=X2S 


(B.21) 


^2   S  (B.22) 


Xo         1      „     ^i 


X2  = 


Xi  I+T2S       I+T2S  (B,23) 


Substituting  equation  B.22  ,    B.23  into  equation  B.21 
we  have : 


X,K^    K-T,X-,S 
Xo  = + 


I+T2S     I+T2S 


(B.24) 
Finally  rearranging  terms  equation  B.24  becomes: 


X3  K^Cl+T-j^S) 

x7'   I+T2S  (B  25) 


92 


D.   CODING  OF  THE  EQUATIONS 
For  Controller  "C" 


Integration- 


YAWE  =  YAWC  -  YAW 

YAWE  =  YAW  -  YAWC 

(YAWE  -  X2) 

DX^  = 

T2 

A  =  K^(X2  +  T^DX2) 

A-X3 
DX3  =  —1: 


X2  =  X2  +  D^2-  ^^^'^ 
X,  =  X3  +  DX3.  DELT 


D  =  X3  +  T^DX3 


For  Controller  "B' 


Integration- 


^AWE  =  YAW  -  YAWC 

YAWE  -  X 
DX2  = _2 

T2 
A  =  K;^(X2  +  T2DX2) 
D  =  X. 


X2  =  X2  +  DX2  •  DELT 
X3  =  X3  +  DX3  •  DELT 


For  Controller  "A' 


Integration- 


YAWE  =  YAW  -  YAWC 
YAWE  -  X_ 


DX, 


X2  =  X2  +  DX2  •  DELT 


D  =  K2(X2  +  Ti    DX2) 
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For  all  of  the  above  cases  the  equations  include  the 
error  detector  and  the  controller  which  are  indicated  in 
Figure  B.4 


YAWC 


YAVJE 


CONTROLLER    'A' 

•B' 
'C 


PLANT 


YAW 


Figure  B.4    General  Scheme  of  Control 
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APPENDIX  C 
RESPONSE  OF  THE  SYSTEM  FOR  REGULAR  SEAS 

//PROGRA  JOB  (????,0356), 'RESEARCH' ,CLASS=A 

//-'MAIN  0RG  =  NPGVM1.  ????P 

//  EXEC  F0RTXCG,PARM.F0RT='0PT(2) ' ,IMSL=DP,REGION=1024K 

//FORT. SYS IN  DD  - 

C 

c 

C  IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS 
C  HAVE  BEEN  OBTAINED  CHANGE  XS(")  TO  X(''-)  AND  DELETE 
c  XU(''-)  ,AND  XL("). 
COMMON  J 
DIMENSION  X(5) 
X(l)=l. 8287125 
X(2)=l. 1652012 
X(3)=10.5659571 
X(4)=ll. 7124157 
X(5)=20. 5683947 
C   CALL  PLANT (X) 

C   IF  ONLY  SIMULATION  IS  WANTED 
CALL  PLANT (X) 
WRITE  (6,25) 
25     FORMAT (IX,'  OPTIMAL  GAINS ',/ ) 

DO  30  1=1,5 
30     WRITE(6,40)I,X(I) 
40     F0RMAT(1X, 'X(' ,12, ' )=' ,F14.7) 
WRITE (6, 50)  J 
50     FORMAT (IX, 'J  =  ',E15.10) 
STOP 
END 

SUBROUTINE  PLANT (XX) 
C   SUBROUTINE  PLANT (XX)  SIMULATES  THE  SHIP 
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COMMON  TDIFF 

REAL" 8  L,L2,L3,L4,L5,L6 

REAL" 8  X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 

REAL-8  TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 

REAL" 8  YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 

REAL" 8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 

REAL" 8  RHO , IZ , FX , FY , MZ , XP , MAS S , DELT , MZ I , RXI , WA , WE 

REAL" 8  DYAWE , YAWE , YAWC , ISE , ISR , LAMDA , D , RYR , RYI , MZR 

REAL-8  Kl , T 1 , T2 , D , X2 , DX2 , S , RX , RY , RZ , TX , TY , TZ , RXR 

REAL-8  T3,T4,X3,DX3,X4 

DIMENSION  XX(5) 
C 

C   CLOSE  LOOP  ANALYSIS  WITH  FILTER 
C 

C   INITIAL  CONDITIONS  FOR  INTEGRATION 
C   SIMULATION  END  TIME  IN  SECONDS 

ETIME=600. 

TIME=0.0 

ICOUNT=1.0 
C   INITIALIZE  THE  COST  FUNCTION 

ISE=0.0 

ISR=0.0 

TDIFF=0.0 

LAMDA=8.128 
C   GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 

K1=XX(1) 

T1=XX(2) 

T2=XX(3) 

T3=XX(4) 

T4=XX(5) 
C      WRITE(6,1010)  K1,T1,T2 

CIOIO   F0RMAT(1X, 'Kl  =',F15.7,'T1  =',F15.7,'T2  =',F15.7) 
C   X, XDOT, Y, YDOT  ARE  FIX  COORDINATES  ON  EARTH 

X=0.0 

Y=0.0 
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XDOT=0.0 

YDOT=0.0 
C   U,UDOT,V,VDOT  ARE  FIX  COORDINATES  ON  SHIP 

V=0.0 

UDOT^O.O 

VDOT=0.0 

YAW=0.0 

R=0.0 

RDOT=0.0 
C   ORDERED  SPEED  IN   FEET/ SEC 
C    38.82  FT/SEC=23  KNOTS 

UC=38.82 
C   AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

U=UC 
C   D  =  RUDDER  ANGLE 

D=0.0 

L=880.5 

L2  =  L^'^^*^2 

L3=L"L"L 

L4=L"L3 

L5  =  L--L4 

L6=L"L5 
C   SEA  DISTURBANCE 

C   FORCES  IN  X,Y  DIRECTION  COMPUTED  IN  FORCES 
C   MOMENTS  IN  Z 

FX=0. 

FY=0. 

MZ  =  0. 
C     RXR=-0.I5744D+05 
C     RXI=-0.19950D+06 
C      RYR=0.52365D+04 
C      RYI=0. 18699D+06 
C     MZR=-0.29870D+08 
C      MZI=-0.35751D+07 
C      RXR=-0.50230D+04 
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C     RXI=0.12712D+05 
C     RYR=0.35290D+04 
C     RYI=-0.31909D+05 
C     MZR=0.38826D+07 
C     MZI=-0.64313D+07 

RXR=0.28540D+04 

RXI=-0.99574D+04 

RYR=-0.85441D+04 

RYI=0.39595D+05 

MZR=-0.13014D+08 

MZI=0.11348D+08 
C     RXR=-0.75642D+04 
C     RXI=0.83497D+04 
C     RYR=0.23379D+05 
C      RYI=-0.81502D+05 
C     MZR=0.28622D+07 
C     MZI=-0.19388D+08 
C      RXR=-0.37916D+04 
C     RXI=0.16381D+04 
C      RYR=-0.76647D+05 
C      RYI=-0.37685D+05 
C     MZR=-0.83915D+07 
C     MZI=-0.53176D+07 

RX  =  DSQRT  (RXR- "2  +  RXI''"2  ) 

RY=DSQRT  (RYR- "2  +  RYI" -'2  ) 

RZ  =  DSQRT  (MZR'-"2+MZI''-"2  ) 

TX  =  DAT AN ( RX I / RXR ) 

T Y = DAT AN ( R Y I / RYR ) 

TZ=DATAN(MZI/MZR) 
C   SIGNIFICANT  WAVE  HEIGHT;  SEA  STATE  1-5,2-10,3-15, 
C   4-17.5,5-22.5  6-27,7-35,8-42,9-60 

WA=17.5 
C   ENCOUNTER  FREQUENCY  .1,. 2, .3, .4, .5, .6, .75, 1.0, 1.5, 2. 5 

WE=0.6 
C   HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  HERE 
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c   AS  PARAMETERS 

RHO=1.9876 

MAS S  =  (  .  0044 )  '*'  ( .  5 '-RH0"L3  ) 

IZ=(0.00028)"(.5"RHO"L5) 

YAWE=0.0 

X2=0.0 

DX2=0.0 

X3=0.0 

DX3=0.0 

X4=0.0 
200  CONTINUE 

S  =  DSQRT  (u*'-'''2+V'''"2 ) 
C   INPUT  YAW  COMMAND 

YAWC=0.0 

IF  (TIME. GE. 0.0)  YAWC=0.0 
C  ERROR  SIGNAL  TO  DRIVE  RUDDER (YAW  ACTUAL -YAW  ORDERED) 
C   (  COMPENSATOR  FILTER  ) 

YAWE=YAW  -  YAWC 

DX2=(YAWE-X2)/T2 

X4  =  Kl-'-(Tl"DX2  +  X2) 

DX3=(X4-X3)/T4 

D=(T3"DX3+X4) 
C   AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 
C   XUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED 
C   FOR  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER 
C   FREQUENCY 

XUDOT=(- .0001)"(  .5"RHO-'L3) 

XU=  (  -  0  .  025 3  )  "  (  .  5 '•'RH0-L2 " S ) 

XUU=  (-0.0003)"  (  .5"RHO-'L2) 

XVR= ( 0 . 00  3  9 ) " ( . 5 "RHO "L3 ) 

XVV= ( - . 00 12 ) " ( . 5 "RHO "L2 ) 

XDD= ( -  0 . 0005 ) " ( . 5 "RH0-L2 " S " "2 ) 
C   LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 
C      YV= ( -  0 . 00  7  5  8 ) " ( . 5  - RHO "L2 - S ) 

YR=(0.0023)"( .5-RHO-L3-S) 
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YD=(0.00145)''(.5"RHO"L2"S^'"'-2) 

YVVR= ( 0 . 0 1 ) " ( . 5 "RHO "L3 / S ) 

YVRR=  (  -  0  .  008  )  -  (  .  5  "RH0'-L4/  S  ) 

YVVV= ( -  0 . 0  3 ) " ( . 5 " RHO "L2 / S ) 

YRRR=  ( 0  .  0 0  3  )  •'  (  .  5  - RHO -'LS  /  S  ) 

YDDD= ( - 0 . 0005 ) " ( . 5 "RH0-L2 - S - -2 ) 
C   YUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE 
C   CHANGED  FOR  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  , 
C   ENCOUNTER  FREQUENCY 
C     YVDOT= ( -0 . 0039 ) " ( . 5 "RH0-L3 ) 

C   SPEED=23  KNOTS,  ENCOUNTER  ANGLE  =  30  ,  ENCOUNTER 
c   FREQUENCY  =0.4 
C     YVDOT=-0.30908D+07 
C     YV=-0.81271D+04 
C      YVDOT=-0.36185D+07 
C     YV=-0.24757D+06 

YVDOT=-0.32890D+07 

YV=-0.11775D+07 
C      YVDOT=-0.23038D+07 
C     YV=-0.18267D+07 
C      YVDOT=-0.59800D+06 
C      YV=-0.13260D+07 
C   MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 

NV=  ( -0  .  00213  )"(  .  5''-RHO"L3 "S  ) 
C      NR= ( -  0 . 00 105 ) " ( . 5 "RHO "L4" S ) 

ND= ( -  0 . 0007 )" ( . 5 "RH0-L3 "S " "2 ) 

NVVR= ( - 0 . 0 15 ) " ( . 5 "RHO "L4 / S ) 

NVRR=  (  -  0  .  008  )  "  (  .  5  "RH0'-L5  /  S  ) 

NVVV= ( 0 . 0 1 ) " ( . 5 "RHO "L3 / S ) 

NRRR=  (  -  0  .  00 6  )  •'-  (  .  5  "RHO "L6  /  S  ) 

NDDD=  ( 0  .  000 1 )  "  (  .  5 '•-RH0"L3  " S " "2  ) 
C   NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE  CHANGED 
C   FOR  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER 
C   FREQUENCY 
C      NRDOT  = ( - 0 . 0  0  0  2  7 ) " ( . 5 -RHO " L5 ) 
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C   SPEED=23  KNOTS,  ENCOUNTER  ANGLE  =  30  ,  ENCOUNTER 

C   FREQUENCY  =0.4 

C  NRDOT=-0.26251D+12 

C  NR=-0.53637D+09 

C     NRDOT=-0.20125D+12 

C  NR=-0.94970D+10 

NRDOT=-0.18671D+12 

NR=-0.46860D+11 
C     NRDOT=-0.14518D+12 
C     NR=-0.87538D+11 
C     NRDOT=-0.37261D+11 
C     NR=-0.69856D+11 
C   REGULAR  WAVE  SEA  STATE 

FX=WA"RX"DCOS(WE"TIME+TX) 

FY=WA"RY'''DCOS(WE'''TIME  +  TY) 

MZ=WA"RZ"DCOS (WE-TIME+TZ ) 
C   U  ACTUAL  SPEED 
C   UC  COMMANDED  SPEED 
C   XP  =  PROPELLER  THRUST 

XP=-XUU'-UC-'"2 
C   EQUATIONS  OF  MOTION 

C      UDOT=(  (XVR  +  MASS)"V'-R  +  XUU"U-'"'-2  +  XVV"V''"*-2 
C     1  +  XDD"D''D  +   FX  +  XP  )  /  (MASS-XUDOT) 

VDOT=(YV"V  +  (YR-MASS"U)"R  +  YD-D 
1  +  YVVR"V""2"R  +  YVRR"V"R""2 
1   +  YVVV-V--3  +  YRRR"R""3  +  YDDD"D""3 
1   +  FY  )/(MASS-YVDOT) 

RDOT=(  NV-V  +  NR'-R  +  ND-D  +  NVVR"V""2"R 
1  +  NVRR"V"R"'-2 

1  +  NVVV-V--3  +  NRRR-R--3  +  NDDD-D-'O 
1  +  MZ  )/(iz-NRDOT) 
C   WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.  2)  GO  TO  50 

GO  TO  300 
C   CONVERT  RADIANS  TO  DEGREES 
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50   YAWDEG=  YAW-57.296 

RDEG  =  R-'-57.296 

RDDEG  =  RDOT'''57.296 

DDEG=D"57.296 

YAWC=YAWC"57.296 

WRITE  (6,100)  TIME,YAWDEG 
C     1   ,UC,U,UDOT,V,VDOT,YAWC,YAWDEG,RDEG,RDDEG,DDEG 

100  F0RMAT(1X,F12.8,1X,F12.8) 
C     1'  FT   XDOT=' ,F8.4, '  FT/ SEC   Y=',F8.2,'  FT 
C     1   YDOT=' ,F8.4, '  FT/ SEC 

C     l',/,2X,'  UC=',F8.4,'  FT/SEC    U=',F8.4,' 
C     1  FT/SEC     UDOT= ' ,F10.6, 
C     1  '  FT/SEC''"2    V=',F8.4,'  FT/ SEC 
C     1  VDOT=' ,F10.6, '  FT/ SEC" -2' 
C     1  ,/,2X, 'YAWC=' ,F8.4, '  DEG   YAW= ' 
C     1    ,F15.7,'  DEG    YAW  RATE=' ,F15.7, ' 
C     1  DEG/ SEC     YAW  ACCEL= ' 
C     1  ,F15.7,'  DEG/SEC--2'  ,/ ,2X,  ' 
C     1  RUDDER  = '  , F15 . 7 , '  DEG  '  ) 
C      WRITE  (6,101)  TIME,DDEG 
ClOl   F0RMAT(1X,F12.8,1X,F12.8) 

IC0UNT=1 
C   TEST  IF  WANT  TO  STOP 

300   IF  (TIME.GE.ETIME)  GO  TO  400 
C   INTEGRATION  STEP  SIZE  DELT 

DELT=1.0 
C   INTEGRATION 

U=U+UDOT'-DELT 

V=V+VD0T'-DELT 

R=R+RD0T"DELT 

YAW=YAW+R'''DELT 

X2=X2+DX2"DELT 

X3=X3+DX3"DELT 
C   CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 

XD0T  =  U"DC0S  (YAW) -V''-DSIN(YAW) 
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YDOT=U"DSIN(YAW) +V"DCOS (YAW) 
X=X+XDOT"DELT 
Y  =  Y  +  YDOT'>DELT 
TIME=TIME+DELT 
IC0UNT=IC0UNT+1 
ISE  =  ISE  +  LAMDA''-YAWE''"'-2 
ISR=ISR  +  D--2 
GO  TO  200 
C   J=TDIFF=  COST  FUNCTION 
400   TDIFF=ISE+ISR 

WRITE (6, 500)  ISE,ISR,TDIFF,K1,T1,T2,T3,T4 
500   FORMATC  '  ,  IX,  '  ISE=  '  ,F15  .  7  ,  '   ISR=  '  ,  F15  .  7  ,  ' 
1  T0TAL=' ,F15.7,2x, 

I  'Kl=' ,F15.7,2X, 'Tl=' ,F15.7,2X, 'T2=' ,F15.7,2X, 
1  'T3=' ,F15.7,2x, 'T4=' ,F15.7) 
RETURN 
END 
//GO. SYS IN  DD  " 
/* 
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APPENDIX  p 
DETERMINATION  OF  OPTIMAL  CONTROLLER  PARAMETERS  FOR  IRREGULAR 

SEAS 

//TRIALl  JOB  (1707,0356), 'RESEARCH' , CLASSIC 
//"MAIN  0RG=NPGVM1.1707P 

//  EXEC  F0RTXCG,PARM.F0RT='0PT(2) ' ,IMSL=DP,REGION=1024K 
//FORT. SYS IN  DD  - 

C   THIS  PROGRAM  WILL  OBTAIN  THE  CONTROLLER  OPTIMAL 
GAINS.  IT  IS  REFERENCED  IN  CHAPTER  5. 
IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS 
HAVE  BEEN  OBTAINED  CHANGE  XS('-)  TO  X(")  AND 
DELETE  XU("),AND  XL(" ) • 

DIMENSION  XS(5),XU(5) ,XL(5) 
XS(1)=0. 655751 
XS(2)=80.5483 
XS(3)=10. 74847 
XS(4)=12.9 
XS(5)=45.09 
XS(I)  IS  THE  STARTING  GUESS 

XL (I)  IS  THE  LOWER  LIMIT  FOR  THE  I ' TH  VARIABLE 
XU(I)  IS  THE  UPPER  LIMIT  FOR  THE  I ' TH  VARIABLE 


xL(i; 

)  =  0.1 

xu(i; 

)  =  2.5 

XL(2; 

)=40.0 

xu(2: 

)=100.0 

XL(3: 

)  =  0.1 

xu(3; 

)=20.0 

XL(4: 

)  =  5.0 

xu(4; 

)=80.0 

XL(5; 

)=60.0 

xu(5: 

)=150.0 

C   A  DESCRIPTION  OF  THE  FOLLOWING  PARAMETERS 
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C   IS  DISCUSSED  IN  BOXPLX 
R=9./13. 
NTA=1000 
NPR=100 
NAV=0 
NV=5 
IP=0 
C   THE  FOLLOWING  STATEMENT  MUST  BE  CHANGED  TO 
C   CALL  PLANT (X) 
C   IF  ONLY  SIMULATION  IS  WANTED 

CALL  BOXPLX (NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , lER ) 
WRITE  (6,25) 
25      FORMAT (IX,'  OPTIMAL  GAINS ',/ ) 

DO  30  1=1,5 
30     WRITE(6,40)I,XS(I) 
40     F0RMAT(1X, 'X(' ,12, ' )=' ,F14.7) 
STOP 
END 

SUBROUTINE  PLANT (XX) 
C   SUBROUTINE  PLANT (XX)  SIMULATES  THE  SHIP 
COMMON  TDIFF 
REAL-8  L,L2,L3,L4,L5,L6 

REAL" 8  X,XDOT,Y,YDOT,U,UDOT,V,VDOT,YAW,R,RDOT 
REAL-8  TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 
REAL" 8  YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 
REAL" 8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL" 8  RHO , IZ , FX , FY , MZ , XP , MASS , DELT 
REAL "8  DYAWE , YAWE , YAWC , I SE , I SR , LAMDA , D 
REAL-8  Kl , Tl , T2 , T3 , T4 , D , X2 , DX2 , X3 , DX3 , X4 , CH ( 11 ) , S 
DIMENSION  XX(5) 
C 

C   CLOSE  LOOP  ANALYSIS  WITH  FILTER 
C 

C   INITIAL  CONDITIONS  FOR  INTEGRATION 
C   SIMULATION  END  TIME  IN  SECONDS 
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ETIME=600. 

TIME^O.O 

IC0UNT=1 
C   INITIALIZE  THE  COST  FUNCTION 

ISE=0.0 

ISR=0.0 

TDIFF=0.0 

LAMDA=8.128 
C   GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 

K1=XX(1) 

T1=XX(2) 

T2=XX(3) 

T3=XX(4) 

T4=XX(5) 
C     WRITE(6,1010)  K1,T1,T2 

CIOIO   F0RMAT(1X, 'Kl  =',F15.7,'  Tl  =',F15.7,' 
C     1   T2=',F15.7   ) 
C   X,XDOT,Y,YDOT  ARE  FIX  COORDINATES  ON  EARTH 

X  =  0.0 

Y=0.0  • 

XDOT=0.0 

YDOT=0.0 
C   U,UDOT,V,VDOT  ARE  FIX  COORDINATES  ON  SHIP 

V=0.0 

UDOT=0.0 

VDOT=0.0 

YAW=0.0 

R=0.0 

RDOT=0.0 
C   ORDERED  SPEED  IN   FEET/ SEC 
C    38.82  FT/SEC=23  KNOTS 

UC=38.82 
C   AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

U=UC 
C   D  =  RUDDER  ANGLE 
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D=0.0 

L=880.5 

L2  =  L— 2 

L3=L"L"L 

L4=L"L3 

L5=L"L4 

L6=L"L5 
C   SEA  DISTURBANCE 

C   FORCES  IN  X,Y  DIRECTION  COMPUTED  IN  FORCES 
C   MOMENTS  IN  Z 

FX  =  0. 

FY=0. 

MZ  =  0. 
C  I SEA  IS  A  SWITCH; I SEA=0  (CALM  WATER)  ISEA=1  (SEA  STATE) 

ISEA=1 
C   HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  HERE  AS 
C   PARAMETERS 

RHO=1.9876 

MAS S  =  ( .  0 044 )  "  (  .  5  -RHO '''LS  ) 

IZ=(0.00028)^M.3"RHO^"L5) 

YAWE=0.0 

X2=0.0 

DX2=0.0 

X3=0.0 

DX3=0.0 

X4=0.0 
200  CONTINUE 

S  =  DSQRT  (U'-"2+V"  "2  ) 
C   INPUT  YAW  COMMAND 

YAWC=0.0 

IF  (TIME. GE. 0.0)  YAWC=0.0 
C   ERROR  SIGNAL  TO  DRIVE  RUDDER(YAW  ACTUAL-YAW  ORDERED) 
C   (  CONTROLLER   FILTER  ) 

YAWE=YAW  -  YAWC 

DX2=(YAWE-X2)/T2 
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X4=K1"(T1"DX2+X2) 

DX3=(X4-X3)/T4 

D=(T3"DX3+X4) 
C   AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 
C   XUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED 
C   FOR  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER 
C   FREQUENCY 

XUDOT= ( - . 0001)" ( . 5-RHO-L3 ) 

XU= ( -  0 . 025  3 ) " ( . 5 "RHO "L2 " S ) 

XUU=  ( -0  .  0003  )'-(  .  5-RHO-L2 ) 

XVR=  (0  .  0039  )"(  .  5"RH0'-L3  ) 

XW=  (  -  .  00 12 )  "  ( .  5  "RHO ''L2 ) 

XDD=  (-0.0005)^^(.  5^'^RH0"L2  "  S  "  "2  ) 
C   LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 

YV= ( -  0 . 00758 ) " ( . 5 "RH0-L2 " S ) 

YR=  (0  .  0023 )"  ( .  5 "RH0^'^L3 "S ) 

YD=  (0  .  00145  )'- ( .  5 "RH0"L2"S--'*-2 ) 

YVVR= ( 0 . 0 1 ) " ( . 5 "RHO "L3 / S ) 

YVRR=  (  -  0  .  008  )  "  (  .  5  "RH0''-L4/  S  ) 

YVVV= ( - 0 . 0  3 ) " ( . 5 "RHO "L2 / S ) 

YRRR= ( 0 . 0  0  3 ) " ( . 5  - RHO - L5 / S ) 

YDDD=  (  -  0  .  0005  )  "  (  .  5  "RH0"L2"S - -'2 ) 
C   YUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED 
C   FOR  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER 
C   FREQUENCY 

C      YVD0T= ( -  0 . 0039 ) " ( . 5 "RH0-L3 ) 
C  SPEED=23  KNOTS,  ENCOUNTER  ANGLE  = 
C  ENCOUNTER  FREQUENCY  =.75 

YVDOT=-2304300.0 
C   MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 

NV= ( - 0 . 002 13 ) " ( . 5 "RH0-L3 " S ) 

NR=  (  -  0  .  00 105  )  ^M  •  5  "RH0^^L4'''S ) 

ND= ( -  0 . 0007 )" ( . 5 "RH0-L3 "S" "2 ) 

NVVR=  (  -  0  .  0 15  )  "  (  .  5  " RHO '-L4  /  S  ) 

NVRR= ( -  0 . 008 ) " ( . 5 "RHO "L5 / S ) 
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NVVV= ( 0 . 0 1 ) " ( • 5 "RHO "L3 / S ) 

NRRR=  (  -  0  .  00 6  )  "  (  .  5  -RHO ''L6  /  S ) 

NDDD=  ( 0  .  000 1 )  "  (  .  5 ''-RHO-'-LS '-S ''"2  ) 
C   NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE 
C   CHANGED  FOR  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  , 
C   ENCOUNTER  FREQUENCY 
C      NRDOT= ( -  0 . 00027 )" ( . 5 »RH0"L5 ) 
C   SPEED=23  KNOTS,  ENCOUNTER  ANGLE  = 
C   ENCOUNTER  FREQUENCY  =.75 

NRDOT=-1.4518E+ll 
C   SETS  SEA  STATE  TO  ZERO 

IF  (ISEA.EQ.l)  GO  TO  30 

FX=0. 

FY=0. 

MZ  =  0. 

GO  TO  35 
C   UNIT  12  HAS  THE  SEA  STATE  DATA  NAMED  CH 
C   IT  MUST  BE  SYNCHRONIZED  BY  TIME 
30    READ  (12)  CH 

FX=CH(3) 

FY=CH(4) 

MZ=CH(8) 
35    CONTINUE 
C   U  ACTUAL  SPEED 
C   UC  COMMANDED  SPEED 
C   XP  =  PROPELLER  THRUST 

XP=-XUU''-UC--2 
C   EQUATIONS  OF  MOTION 

C      UDOT=(  (XVR  +  MASS)"V-''R  +  XUU"U"'-2  +  XVV"V"''2 
C     1  +  XDD'*'D"D  +   FX  +  XP  )  /  (MASS-XUDOT) 

VDOT=(YV"V  +  (YR-MASS"U)"R  +  YD-D  +  YVVR''-V"-'2"R 
1   +  YVRR"V"R""2 

1   +  YVVV"V""3  +  YRRR-R--3  +  YDDD"D""3 
1   +  FY  )/(MASS-YVDOT) 

RDOT=(  NV"V  +  NR-R  +  ND-D  +  NVVR"V""2"R 
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1   +  NVRR"V"R''"2 

1  +  NVVV"V'"'-3  +  NRRR"R""3  +  NDDD"D""3 
1  +  MZ  )/(IZ-NRDOT) 
C   WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.il)  GO  TO  50 
GO  TO  300 
C   CONVERT  RADIANS  TO  DEGREES 
50   YAWDEG=  YAW''-57.296 
RDEG  =  R''-57.296 
RDDEG  =  RDOT'-57.296 
DDEG=D"57.296 
YAWC=YAWC'-57.296 
C      WRITE  (6,100)  TIME,XP,X,XDOT,Y,YDOT 
C     1   ,UC,U,UDOT,V,VDOT,YAWC,YAWDEG,RDEG,RDDEG,DDEG 
100  F0RMAT(1X, 'TIME=' ,F8.3, '  SEC   XP=',F10.2,'  LBF 
1   X=',F8.2, 

1'  FT   XDOT=' ,F8.4, '  FT/ SEC   Y=',F8.2,'  FT   YDOT= ' 
1  ,F8.4,'  FT/SEC 

l',/,2X,'  UC=',F8.4,'  FT/SEC    U=',F8.4,'  FT/SEC 
1  UDOT=' ,F10.6, 

1  '  FT/ SEC" "2    V=',F8.4,'  FT/ SEC    VDOT= ' , FIO . 6 , 
1  •  FT/SEC"'-2'  ,/,2X,  •YAWC=' ,F8.4,  '  DEG   YAW=  ' 
1    ,F15.7,'  DEG    YAW  RATE=' ,F15.7, '  DEG/ SEC 
1   YAW  ACCEL=' 

1  ,F15.7,'  DEG/ SEC'"'-2' ,/ ,2X, 'RUDDER  ='  ,F15.  7,  ' 
1  DEG  • , / ) 
I COUNT =1 
C   TEST  IF  WANT  TO  STOP 

300   IF  (TIME.GE.ETIME)  GO  TO  400 
C   INTEGRATION  STEP  SIZE  DELT 

DELT=1.0 
C   INTEGRATION 

U=U+UDOT"DELT 
V=V+VDOT-'-DELT 
R=R+RDOT"DELT 
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YAW=YAW+R"DELT 

X2  =  X2  +  DX2'''DELT 

X3  =  X3+DX3"-DELT 
C   CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 
C      XDOT  =  U"DCOS(YAW)-V''DSIN(YAW) 
C      YDOT  =  U''DSIN(YAW)+V"DCOS(YAW) 
C      X  =  X+XDOT'-DELT 
C      Y=Y+YDOT'-DELT 

TIME=TIME+DELT 

IC0UNT=IC0UNT+1 

ISE=ISE  +  LAMDA"YAWE-'"'-2 

ISR=ISR  +  D''""2 

GO  TO  200 
C   J=TDIFF=  COST  FUNCTION 
400   TDIFF=ISE+ISR 

WRITE (6, 500)  TDIFF,K1,T1,T2,T3,T4 
500   FORMAT ('  ',1X,'TDIFF  =',F15.7,'   Kl  =',F15.7,' 
1   Tl  =' ,F15.7,2X, 
1  'T2=' ,F15.7,2X, 'T3=' ,F15.7.2X, 'T4=' ,F15.7) 

REWIND  12 

RETURN 

END 
C      BETWEEN  LINE  249 (END)  AND  THE  FOLLOWING  LINE 
C      (//GO.  SYS  IN  DD  -'OWE  HAVE  TO  INCLUDE  BOXPLX. 
//GO. SYS IN  DD  " 

//GO.FT12F001  DD  DISP=SHR,DSN=MSS . S2160 . A213 
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APPENDIX  E 
RESPONSE  OF  THE  SYSTEM  FOR  IRREGULAR  SEAS 

//PROGRA  JOB  (????,0356) , 'RESEARCH' ,CLASS=B 
//'•-MAIN  0RG  =  NPGVM1.????P 
//  EXEC  FRTXCLGP,IMSL=DP,REGION=1024K 
//FORT. SYS IN  DD  - 

C   IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS 
C   HAVE  BEEN  OBTAINED. 
DIMENSION  XX(5) 
C   OPTIMAL  GAINS  FOR  CONTROLLER 
XX(1)=2. 45967680 
XX(2)=88. 2797241 
XX(3)=50. 5678864 
XX(4)=5. 27039050 
XX(5)=95. 3189392 
C   THE  SUBROUTINE  PLANT  SIMULATES  THE  SL-7  CONTAINERSHIP 

CALL  PLANT (XX) 

WRITE (6, 25) 
25     FORMAT (IX, 'OPTIMAL  GAINS' ,/) 

DO  30  1=1,5 
30     WRITE (6, 40) I, XX (I) 
40     F0RMAT(1X, 'XX(' ,12, ' )=' ,F14.7) 

STOP 

END 
C 
C   SUBROUTINE  PLANT (XX)   SIMULATES  THE  SHIP 

SUBROUTINE  PLANT (XX) 

COMMON  TDIFF 

REAL" 8  L,L2,L3,L4,L5,L6 

REAL '-8  X ,  XDOT  ,  Y ,  YDOT  ,  U ,  UDOT  ,  V ,  VDOT  ,  YAW ,  R ,  RDOT 

REAL" 8  TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 

REAL" 8  YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 
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REAL" 8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 

REAL "8  RHO , I Z , FX , FY , MZ , XP , MAS  S , BELT 

REAL "8  DYAWE , YAWE , YAWC , I SE , I SR , LAMDA , D 

REAL"8  K1,T1,T2,D,X2,DX2,S,CH(11) ,DX3,X3,X4 

DIMENSION  XX (5) 
C 

C   CLOSE  LOOP  ANALYSIS  WITH  FILTER 
C 

C   INITIAL  CONDITIONS  FOR  INTEGRATION 
C   SIMULATION  END  TIME  IN  SECONDS 

ETIME=600. 

TIME=0.0 

I COUNT =1 
C   INITIALIZE  THE  COST  FUNCTION 

ISE=0.0 

ISR=0.0 

TDIFF=0.0 

LAMDA=4.2 
C   GAIN  COEFFICIENTS 

K1=XX(1) 

T1=XX(2) 

T2=XX(3) 

T3=XX(4) 

T4=XX(5) 
C   X,XDOT,Y,YDOT  ARE  FIX  COORDINATES  ON  EARTH 

X=0.0 

Y  =  0.0 

XDOT=0.0 

YDOT=0.0 
C   U,UDOT,V,VDOT  ARE  FIX  COORDINATES  ON  SHIP 

V=0.0 

UDOT=0.0 

VDOT=0.0 

YAW=0.0 

R=0.0 
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RDOT^O.O 

YAW=0.0 
C   ORDERED  SPEED  IN   FEET/ SEC 
C    54.01  FT/SEC=32  KNOTS 

UC=38.81 
C   AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

U=UC 
C   D  =  RUDDER  ANGLE 

D=0.0 

L=880.5 

L2=L""2 

L3=L"L"L 

L4=L"L3 

L5  =  L'-L4 

L6=L"L5 
C   SEA  DISTURBANCE 

C   FORCES  IN  X,Y  DIRECTION  COMPUTED  IN  FORCES 
C   MOMENTS  IN  Z 

FX=0. 

FY=0. 

MZ  =  0. 
C   ISEA  IS  A  SWITCH;  ISEA=0(CAL  WATER) ISEA=1 (SEA  STATE) 

ISEA=1 
C   HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  HERE  AS 
C   PARAMETERS . 

RHO=1.9876 

MASS= ( . 0044) " ( . 5-RHO-L3 ) 

IZ= (0 . 00028 )"( " 5-RHO-L5 ) 

YAWE=0.0 

X2=0.0 

DX2  =  0.0 

X3  =  0.0 

DX3=0.0 

X4=0.0 
200  CONTINUE 
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S  =  DSQRT(U--2  +  V''"2) 
C   INPUT  YAW  COMMAND 

YAWC=0.0 

IF  (TIME. GE. 0.0)  YAWC=0.0 
C   ERROR  SIGNAL  TO  DRIVE  RUDDER (YAW  ACTUAL- YAW  ORDERED) 
C   (  COMPENSATOR  FILTER  ) 

YAWE=YAW  -  YAWC 

DX2=(YAWE-X2)/T2 

X4=K1"(T1"DX2+X2) 

DX3=(X4-X3)/T4 

D=(T3''DX3  +  X4) 
C   AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 
C 

C   XUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED 
C   FOR  DIFFERENT  ENCOUNTER  ANGLE  AND  SPEED. 
C      XUDOT  = ( - . 0  0  0 1 ) " ( . 5 "RHO "L3 ) 

XUU= ( -  0 . 0003 )" ( . 5 "RH0-L2 ) 

XVR= ( 0 . 00  3  9 ) " ( . 5 "RHO "L3 ) 

XVV=  ( -  .  00 12 )  "  (  .  5  "RH0^'^L2 ) 

XDD=  (  -  0  .  0005  )  "  (  •  5 '•-RH0-L2 " S "  "2 ) 
C   LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 

YV=  (  -0  .  00758  )^M  .  5"RH0^''L2"S ) 

YR= (0 . 0023 )" ( . 5 "RH0-L3 "S ) 

YD=  ( 0  .  00 145  )  "  (  .  5  "RHO "L2 '-S "  "2  ) 

YVVR= ( 0 . 0 1 ) " ( . 5 "RHO "L3 / S ) 

YVRR= ( -  0 . 008 ) " ( . 5 "RH0-L4/ S ) 

YVVV=  (  -  0  .  0 3  )  "  (  .  5  "RHO "L2  /  S  ) 

YRRR=  ( 0  .  0 0 3  )  "  (  .  5  - RHO '•L5  /  S  ) 

YDDD=  (  -  0  .  0005  )  "  (  .  5 '•-RH0-L2 - S -"''2  ) 
C   YUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED 
C   FOR  DIFFERENT  ENCOUNTER  ANGLE  AND  SPEED. 
C 
C      YVDOT=(-0.0039)"( .5-RHO-L3) 

YVDOT=-2304300.00 
C   MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 
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NV=  (  -  0  .  00213  )"  (  .  5  "RH0^*a3  "S ) 

NR=  ( -0  .  00105  )"(.  5'-RHO^''L4^'^S  ) 

ND=  (  -  0  .  0007 )"(  .  5"RHO"L3^'^S^"'''2 ) 

NVVR=  (  -  0  .  0 15  )  "  (  .  5  "RHO '-L4  /  S  ) 

NVRR= ( -  0 . 008 ) " ( . 5 "RHO "L5 / S ) 

NVVV= ( 0 . 0 1 ) " ( . 5 "RHO "L3 / S ) 

NRRR=  (  -  0  .  0 0 6  )  •'  (  .  5  -RHO "L6  /  S  ) 

NDDD= ( 0 . 000 1 ) " ( . 5 "RH0-L3 " S " "2 ) 
C   NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE  CHANGED 
C   FOR  DIFFERENT  ENCOUNTER  ANGLE  AND  SPEED. 

C 

C     NRDOT=  (  -0  .  00027  )"(  .  5'*-RHO''L5  ) 

NRDOT=-1.5096E+ll 
C   SETS  SEA  STATE  TO  ZERO 

IF  (ISEA.EQ.l)  GO  TO  30 
FX=0. 
FY=0. 
MZ  =  0. 
GO  TO  35 
C   UNIT  12  HAS  THE  SEA  STATE  DATA  NAMED  CH 
C   IT  MUST  BE  SYNCHRONIZED  BY  TIME 
30     READ  (12)  CH 
FX=  CH(3) 
FY=  CH(4) 
MZ=  CH(8) 
35     CONTINUE 
C   U  ACTUAL  SPEED 
C   UC  COMMANDED  SPEED 
C   XP  =  PROPELLER  THRUST 

XP=-XUU"UC"'-2 
C   EQUATIONS  OF  MOTION 

C     UDOT=(  (XVR  +  MASS)"V"R  +  XUU"U''"'-2  +  XVV-V--2 
C     1  +  XDD-D-D  +  FX   +  XP  ) / (MASS-XUDOT) 

VDOT=(YV"V  +  (YR-MASS"S)''-R  +  YD-D  +  YVVR"V'"'-2"R 
1   +YVRR"V-''R""2 
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1   +  YVVV'''V"''3  +  YRRR'-R-"'-3  +  YDDD'-D'"'-3 
1   +  FY  )/(MAS-YVDOT) 

RDOT=(  NV-'V  +  NR'-R  +  ND-'-D  +  NVVR"V--'-2''R 
1  +  NVRR"V"R""2 

1  +  NVVV-V--3  +  NRRR"R""3  +  NDDD-D--3 
1  +  MZ  )/(IZ-NRDOT) 
;   WHEN  TO  PRINTOUT 

IF  (IC0UNT.EQ.2  )  GO  TO  50 

GO  TO  300 
;   CONVERT  RADIANS  TO  DEGREES 
50   YAWDEG=  YAW-'57.296 

RDEG  =  R--57.296 

RDDEG  =  RDOT'-57.296 

DDEG=D"57.296 

YAWC=YAWC"57.296 

WRITE  (6,100)  TIME,YAWDEG 
100  F0RMAT(1X,F12.8,1X,F12.8) 

IC0UNT=1 
;   TEST  IF  WANT  TO  STOP 
300   IF  (TIME.GE.ETIME)  GO  TO  400 
;   INTEGRATION  STEP  SIZE  DELT 

DELT=1. 
;   INTEGRATION 

U=U+UDOT'-DELT 

V=V  +  VDOT--DELT 

R=R+RDOT"DELT 

YAW=YAW+R"DELT 

X2=X2+DX2"DELT 

X3  =  X3  +  DX3-''DELT 
:   CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 

XDOT=U"DCOS(YAW)-V"DSIN(YAW) 

YDOT=U'-DSIN(YAW)+V"DCOS(YAW) 

X=X+XDOT"DELT 

Y  =  Y  +  YDOT'-DELT 

TIME=TIME+DELT 
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IC0UNT=IC0UNT+1 
ISE=ISE  +  LAMDA"YAWE'"'-2 
ISR=ISR  +  D""2 
GO  TO  200 
C   J=TDIFF=  COST  FUNCTION 
400   TDIFF=ISE+ISR 

WRITE(6,500)  ISE,ISR,TDIFF 
500   F0RMAT('l' ,5X, 'ISE=' ,F15.7, '   ISR= ' ,F15 . 7 , ' 
1  T0TAL= ' ,F15.7) 
STOP 
END 
//GO. SYS IN  DD  " 

/" 

//GO.FT12F001  DD  DISP=SHR,DSN=MSS . S2160 . A211 
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